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SYNOPSIS 

Singular Perturbation, now is a maturing mathematical subject 
with a fairly long history and a strong promise for continued 
important applications throughout science and engineering. A 
Singular Perturbation problem is well defined as one in which no 
single asymptotic expansion is uniformly valid throughout the 
interval, as the perturbation parameter e — » 0. The singular 
perturbation problems find place in many area of engineering and 
applied mathematics, for instance: fluid mechanics, fluid 
dynamics, quantum mechanics, plasticity, chemical reactor theory 
and many other problems of fluid motion. A few notable examples 
are boundary layer and WKB problems. The problems to be solved by 
means of asymptotic and numerical analysis are becoming 
progressively more and more complicated, j 

The aim of this thesis is to derive some simple and efficient 
numerical methods using splines for solving singular perturbation 
problems which are easy to implement and are not costly in terms 
of computer time also. 

Chapter I of the thesis contains general introduction which 
includes a brief survey of asymptotic and numerical analysis of 
Singular perturbation problems. A summary of recent methods and 
that of the present thesis is presented. 

In chapter II, we first present a numerical method namely 
Method 1, for Linear Singularly Perturbed Initial Value Problems 
of the form: 

ey' + P(x,e)y = q(x,e); y(a) = a, a * x * b, 



xin 


where e is a small positive parameter, 0 < e « 1, ocisa given 
constant and p(x,e) > o and q(x,e) are sufficiently smooth 
functions. The method consists of dividing the inner and outer 
region and solving the inner region problem by a fourth order 
cubic spline after rescaling the inner region. The solution of the 
reduced problem provides the solution of the outer region problem. 
S-inee "bhe point of division of inner and outer region is arbitrary 
and also it is observed that the solution obtained is 
discontinuous at the point of division. In an attempt to overcome 
these drawbacks, a more general method namely Method 2 is derived, 
which is also applicable to Non-Linear Problems of the form: 

ey /= f (x,y,e) ; y(a)= a, a x =s b, 

where is non-zero on (a,b) and has only one stable root 

for all values of x and y. The method consists of a numerical 
cutting point technique. This method also takes care of continuity 
at the cut point. Both these methods provide global fourth order 
approximation to the exact solution and does not depend on 
asymptotic expansion. 

In chapter III, A class of Non-Linear Two-Point Singularly 
Perturbed Boundary Value Problems is considered. 

ey ,7 (x) + [p(y(x) ) ] ' + q (x,y(x)) = r(x) ; y(a) = a, y(b) = /3; 
here, 0 < e « 1; a and 0 are given constants; p(y) , q(x,y) , r(x) 
are assumed to be sufficiently differentiable functions. 
Furthermore, we assume that the problem has a solution which 
displays a boundary layer of width o(e) at x = a for small value 



of £• The original boundary value problem is reduced to an 
asymptotically equivalent following initial value problem by an 
initial value technique: 

ey' (x) + p(y (x) ) + Q(x) = R(x) + K; y(a) = a, K a constant. 

Then it is solved in the inner region by cubic spline method for 
inner region employed in Method 2 developed in chapter 2. The 
solution of the reduced problem: 

[P(Y 0 (x))]' + q(x,y Q (x)) - r(x) ; y Q (b) = /3, 
provides the solution in the outer region. 

In chapter IV, we first consider the following class of general 
Linear Singularly Perturbed Boundary Value Problems without first 
derivative term. 

c y" = q(x) y + r(x) ; y(a) - a Qf y(b) = oc ± , 
where q(x) > 0 and r(x) are sufficiently smooth functions. 

Ahlberg's cubic spline in terms of moments is considered. A second 
order variable mesh finite difference scheme is developed using the 
continuity conditions of second derivatives of the spline at the 
interior nodes. Then the variable mesh scheme is generalised for 
the following class of problems with first derivative term: 

£ y"= p (x) y # + q(x) y + r(x); y(a) = oc Q , y(b) = a ±t 
where p(x) , g(x) and r(x) are sufficiently smooth functions. 

In deriving the difference scheme for above class of 
problems, we have taken some second order approximations of first 
derivative terms in terms of solution of the problem at the nodal 
points. 
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In Chapter V, The following class of general Semi-Linear Singularly 

Perturbed Boundary Value Problems is considered. 

cy" = f(x,y); y(a) = a, y(b) = /3; a < x < b # 0 < c « 1. 

It is known that, under suitable assumptions, one expects the 
solution to behave qualitatively as the solution of the following 
problem: 

ey" = p(x) y + q(x) ; y(a) = a , y(b) = fi; p(x) 2 : p > 0, 
where p(x) and q(x) are sufficiently smooth functions. 

Therefore, first we consider the above linear singularly perturbed 
problem. A family of variable mesh methods based on cubic spline 
approximations which gives third order approximations to the 
solution of the linear problem. Then these methods are used for 
solving semi-linear problems using a quasi-linearisation 
technique. When a parameter namely mesh ratio parameter present in 
the methods is taken to be unity, the family of third order 
methods reduces to a family of fourth order methods, and in 
addition, if other parameter is also taken to be unity, the family 
of methods reduces to a well known Numerov's method for 
corresponding regular problem. » 

In Chapter VI, the variable mesh methods based on cubic spline 
approximations developed in the previous chapter are generalised 
for the general Non-linear Singularly Perturbed Boundary Value 
Problems: 

ey" = f(x,y,y'); y(a) = A, y(b) = B; x e (a,b) , 0 < e « 1. 



CHAPTER i 


Introduction 

1.1 GENERAL : 

Singular Perturbation, now is a maturing mathematical subject 
with a fairly long history and a strong promise for continued 
important applications throughout science and engineering. We give 
the brief definition of a singularly perturbation problem in its 
simplest and most widely used form. Consider a problem P £ in some 
differential model, depending on a small positive parameter e, 
where 0 < c « 1. Under some conditions, a solution y e (x) of the 
problem P £ can be constructed by the well known method of 
perturbation i.e. as a power series in e with first term y Q being 
the solution of the problem P Q , which is obtained by putting c 
equal to zero in the problem p £ . Under the happiest circumstances, 
this perturbation method leads to altogether satisfactory results. 
This series cannot often be presumed to uniformly converge, 
particularly for small values of e, in the entire interval, when 
such an expansion converges as e — > 0, uniformly in x, one speaks 
of 'Regular Perturbation Problem'. On the other hand, when y £ (x) 
does not have a uniform limit in x as e — » 0, this straight forward 
perturbation method fails and as a consequence of the 
non-uniformity, one may miscalculate or even loss essential 
results, one then speaks of 'Singular Perturbation Problem'. A 
Singular Perturbation problem is well defined as one in which no 
single asymptotic expansion is uniformly valid throughout the 
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interval, as c — > 0. The prototype of Singular Perturbation Problem 
is Prandtl's boundary layer theory, which provides a beautiful 
mathematical tool for the investigation of a practical problem like 
Navier— stokes equations with large reynold number, one of the most 
striking example of singular perturbation. By definition, the 
boundary layer is a narrow region where the solution of a 
differential equation changes rapidly and the thickness of boundary 
layer must approach zero as c tends to zero. 

The singular perturbation problems find place in many area of 
engineering and applied mathematics, for instance: fluid mechanics, 
fluid dynamics, quantum mechanics, plasticity, chemical reactor 
theory, magneto hydrodynamics, reaction-diffusion processes and 
many other problems of fluid motion. A few notable examples are 
boundary layer and WKB problems. The problems to be solved by means 
of asymptotic and numerical analysis are becoming progressively 
more and more complicated. 

The aim of this chapter is to present a brief survey of 
asymptotic and numerical analysis of Singular Perturbation 
Problems. A brief survey on the solution behavior of singular 
perturbation problems is also included. A summary of some recent 
methods is presented. Finally, the summary of the thesis is 
presented. 

1.2 ASYMPTOTIC ANALYSIS OF SINGULAR PERTURBATION PROBLEMS : 

In the first part of the century, although Birkhoff (1908) , 
Langer(1931) and others have done work on asymptotic analysis of 
Linear Ordinary Differential Equations and significant amount of 
work on turning point problems being done by physicists 
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Wentzel (1926) , Krammers (1926) , Brillouin(1926) and others, but 
Friedrichs and Wasow seem to be the first mathematicians to study 
asymptotic solution of singularly perturbed boundary value 
problems. They first use the term 'Singular Perturbation' in print 
in the title of Friedrichs and Wasow(1946). Their work was 
motivated by an analysis of the edge effect for buckled plates. 
Soon afterwards, many authors like Tikhonov ( 1948 ) , Levinson(1950) , 
and Vasileva and Volosov(1967) , began studying related problems. 
Levinson started study of wide range of important topics in 
Singular Perturbation and made contribution to the area together 
with number of young collaborators including Coddington, Davis, 
Flatto, Haber and Levin. The Russian school also did outstanding 
work on many subjects including boundary layer methods — Vishik 
and Lyusternik(1961) . From around 1950, fluid dynamists solved some 
very interesting problems like linoleum-rolling problem — 
Carrier (1953) and low reynold number flow past bodies — 
Kaplaun(1957) . At Caltech's Guggenheim Aeronautical Laboratory, 
Lager strom, Cole, Latta, Van Dyke, Kaplaun and others became 
equally involved in asymptotic expansion procedures for more 
general singular perturbation problems. An over-simplified matching 
procedure is presented in the book of Van Dyke (1964) . The straight 
forward recipe he provided made it easy for tremendous variety of 
scientists to learn the rudiments of matching and to the important 
problem in their own disciplines. The basic idea, much as in 
Friedrichs lectures and Erdeyli (1956) lectures, involved an 
asymptotic matching of the inner and outer expansions at the edge 
of the boundary layer, where they should be appropriate. Cole (19 68) 
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stressed the limit process expansions and two timings m a context 
fa.r broader than fluid mechanics. Indeed the results obtained 
through matching generally coincided with those known through the 
intuitive folkways of the various fields. Wasow(1965) placed the 
singular perturbation in the contexts of the analytic theory of 
differential equations. By 1970, courses in perturbation methods 
became common in science, engineering and applied mathematics 
departments, and inevitably a string of text books and high level 
monograph began to appear. They include, Erdelyi (1956) , Bellman 
and Cooke ( 1963) , Bellman(1964) , Van Dyke (19 64 ) , Kaplaun (1967) , 
Carrier and Pearson(1968) , Ames(1972) ,El'sgol'ts and Norkin (1973 ) , 
Dingle(1973) , Willoughby ( 1974) , O'Malley (1974) , Aziz (1975), Brauner 
et al. (1977), Driver ( 1977 ) , Bender and Orszag(1978) , 
Eckhaus (1973, 1979) , Na(1979), Childs et al.(1979), Hughes (1979) , 

Hemker and Miller(1979) , Axelsson et al. (1979) , Meyer and 
Parter(1980) , Doolan et al.(1980), Nayfey(1981) , Kevorkian and 
Cole(1981), Miranker (1981) , Eckhaus and de Jager(1982), 

Ardema ( 1983 ) , Verhlust (1979 , 1983) , Miller (1975, 1976, 1980, 1982) . 

In O'Malley 's (1982) book review, an outline of history of 
asymptotic methods for singular perturbation problems is given. 

1.3 SOLUTION BEHAVIOR OF SINGULAR PERTURBATION PROBLEMS ; 

Many authors have studied the behavior of singularly perturbed 
problems also. 

Hoppensteadt (1966) has discussed the behavior of the singularly 
perturbed problems on infinite interval. They considered the 
initial value problems of the form 
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x' = f (t,x,y,e) ; x(0) = x 

( 1 . 1 ) 

ey 7 = g(t,x,y,e); y(0) = y Q 

where 0 < e << 1 and x, f are real K-dimensional vectors with 
components x == (x^^,...,^) and f = (f ^ , f , . . . , f^) , respectively 
and y and g are real J-dimensional vectors with components y = 
( y l ,Y 2 / * * • ,y j) an< ^ 9 “ ’ *9j) ' respectively. The purpose of 

this paper was to investigate the behavior of solutions of (1.1) as 
c — » 0 , for t Q 3 t ^ oo. 

Dorr and Parter (1970) have discussed the asymptotic behavior as 
e — » 0 of solutions u(t) = u(t,e) and v(t) = v(t,c) to nonlinear 
boundary value problems of the form 

u" = f(t,u,v) ; 0 < t < 1, u ( 0) = u(l) = 0 (1.2) 

ev" + g(t,u,u')v' - C(t ,u,u') = 0 

(1.3) 

0 < t < 1, v(0) = v Q ; v ( 1 ) = v ± 

where, they have assumed 0 ^ v Q £ and C(t,u,u') £ 0. They were 

particularly interested m problems in which there is exactly one 
interior turning point for the equation (1. 3) , that is, for each 
e > 0 there is unique point a e (0,1) such that g (<x,u (a) ,u' (a) ) =0 
and g(t,u(t) ,u' (t) ) changes sign in a neighbourhood of t = a. In 
order to study the asymptotic behavior of the solutions they 
assume d the following conditions : 

(i) 0 s v o s v i' 

(ii) C(t,u,u') £ 0 

(iii) f(t,u,u'), g(t,u,u') and C(t,u,u / ) are continuous in all 
variables 
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(iv) There exists a continuous function f Q (t,v) such that 
| f(t,u,v) | * f 0 (t , v) for t e [0,1] and v € [0,v 1 ]. 

Hoppensteadt ( 197 1) , investigated the properties of solution of 
ordinary differential equations with a snail parameter. In this 
paper they included the case in which the boundary conditions also 
depend on e i.e. x(t Q ) = a(e) and y(t Q ) = 0(e). They also extended 
these results to the boundary value problems. 

Cohen (197 3) has discussed the existence and asymptotic behavior for 
small £ > 0 of the solution of the nonlinear two-point boundary 
value problem 

ey" + f(x,y,y / )y / =0, 0 < x < l (1.4) 

with 

y' (0) - ay (0) = A 0 ( a > 0 ) (1.5a) 

y' (1) + by ( 1) = B > 0 ( b > 0 ) (1.5b) 

They have imposed the following conditions on the nonlinear f m 
(1.4) 

(i) f(x,y,y / ) is continuously differential in the region 
R = { (x,y, y 7 ) : 0 * x * l, 0* y < B/b, y' * 0} 

(ii) y ± * y 2 and y^ ^ y' 2 imply that 

f (x,y 1 ,y^) ^ f(x,y 2 ,y^) on R 
(m) f(x,y,y') £ 0 > 0 on R 

(iv) There exists a constant k such that for all (x,y,y') € R 
I f (x,y,y') - f(x,z,z'} I a k( |y - z| + I y'-z'l). 

Under these conditions, they proved that there exists a solution 
y(x,e) of (1.4), (1.5a & 1.5b) such that y(x,e) tends to B/b and 
y'(x,e) tends to zero on any subinterval 0 < 5 ^ x =s l. The entire 
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analysis is based on the so-called 'shooting method' for ordinary 
differential equations. 

Kriess and Parter(1974) have discussed the behavior of the 
solutions of the boundary value problems with turning points: 

ey"(x) + f (x,e)y' (x) + g(x,e)y(x) = 0 (1.6) 

for -a =s x * b with y(-a) = A and y(b) = B (1.7) 

where a,b > 0, e > 0 and f(x,e) has a simple zero in [-a,b]. 

Sannuti (1975) has considered a class of two-point boundary value 
problems 

x' = g- (x, t) + B (t) z + C (t)u 

( 1 . 8 ) 

Ay # = g 2 (x,t) + B 2 (t)z + C 2 (t)u 

with x ( 0) = x Q and z(0) = z Q (1.9) 

where x and z are n and m dimensional state vectors, respectively, 
u is an r-dimensional control vector, and A is a nonnegative scalar 
parameter and g^, g 2 , B^ , B 2 , C^, C 2 are assumed to be infinitely 
differentiable in all their arguments m an appropriately defined 
domain, which arise in fixed final time free end point optimal 
control problems. 

Harris (1976) has described the applicability of differential 
inequalities to singular perturbation problems by studying the 
model nonlinear singular perturbation problem 

ey" + yy 7 - y = 0 ; O^x^l (1.10) 

y(0) = A and y(l) = B (1.11) 

whose solution exhibit a wide variety of interesting behavior. 
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Kopell and Parter(1981) have discussed the analysis of the problem 

ey"(t,e) = [y 2 (t,e) - t 2 ]y'(t,c) (1.12) 

with y(-l,e) = A, and y(0,e) = B (1.13) 

based entirely on priori estimates and the 7 shooting method' . 
Grasman and Maktowsky (1977) have employed a variational formulation 
of the problem (1.12) with (1.13) to resolve the question of the 
number of location of the boundary layers as well as to uniquely 
determine the asymptotic expansion of the solution. These results 
are extended to analogous problems for partial differential 
equations with turning points. 

Nipp(1983) has studied an extension of Tikhonov's theorem in 
singular perturbations. They considered the autonomous system of 
ordinary differential equations. 

Finden(1983) has presented an asymptotic approximation for solving 
singular perturbation problems. The given system 


x' = f(t,x,y,e) 
ey' = g(t / x # y / e) 


(1.14) 


together with initial conditions 

x(0) = a, y (0) - /3 (1-15) 

where x, f and a are m-dimensional vectors and y, g, are 
n-dimensional vectors, is replaced by the following system which is 
not stiff numerically. 


x' = f (t,x,y,e) 
y' = h(t,x,y,e) 


for some function h. 


(1.16) 
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They have proved that the solution to (1.16) is a second order 
approximation for the quantity g(x,y,t)/e is given by 

3 — e Y,t ^ = h(x,y,t,e) + 0(e 2 ). (1.17) 

Dadfar et al. (1984) has constructed and analysed a power series 
expansion in the damping parameter e of the limit cycle U(t,e) of 
the free Van der Pol equation 

U" + e(U 2 - l) U' + u = 0 (1.18) 

Grasman ( 1979 ) has discussed similar results for a class of elliptic 
singular perturbation problems. 

Howes (1979) has discussed the boundary layer behavior of singularly 
perturbed initial value problems and the third order boundary value 
problems in Howes (1982). 

Howes (1984) presented some results and sufficient condition in 
order that the solution of the problem (scalar and vector form) 

ey" = F(y)y' + g(x f y) , a s x ^ b (1.19) 

with y(a,e) = a and y(b,e) = ft, (1.20) 

display boundary layer behavior at an end point. 

Axelsson and Carey (1985) considered a class of singularly perturbed 
boundary value problem with various type of boundary conditions and 
examined the boundary layer behavior of the solution. By 
constructing a more regular modified problem with a correction 
term, they discussed how to deal with the boundary layer separately 
and gave error estimates which are uniform m singular perturbation 
parameter. 
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1. 4 NUMERICAL ANALYSIS OF SINGULAR PERTURBATION PROBLEMS : 

The asymptotic analysis treats restrictive class of problems and 
require the problem solver to have some understanding of the 
behavior of the solution expected. As far as an asymptotic analysis 
is valid and a few terms in the asymptotic expansion describe the 
solution sufficiently accurate, one can rely on standard techniques 
to obtain the solutions. As and when the asymptotic analysis is 
difficult to handle or perform badly, then one can usually look for 
an alternative approach and because of this numerical analysis for 
solving singular perturbation problems came into existence. We 
would say that numerical analysis tries to provide quantitative 
information about a particular problem, whereas asymptotic analysis 
tries to gain insight about the qualitative behavior of a family of 
problems and only semi quantitative information of any particular 
member of the family of problems. Numerical methods are intended 
for broad class of problems and are intended to minimize demand 
upon the problem solver. Pearson (1968a) was perhaps the first to 
attempt to give numerical solution of singular perturbation 
problems. After him, in later years, many efficient numerical 
methods were developed, which can be broadly classified as finite 
difference methods, finite element methods and spline approximation 
methods . 

1.4.1 FINITE DIFFERENCE METHODS : 

Pearson (1968a) was perhaps the first to attempt something like 
net adjustments in the difference schemes while treating singular 
perturbation problems. Basically, his idea was to use variable mesh 
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width. Based on the classical three-point difference scheme for a 
non-uniform mesh, Pearson gave the numerical solution of a great 
variety of singular perturbation problems of the form 

L c [y] = e y"(x) + a (x) y' (x) + b(x)y(x) = f(x), (1.21) 
for -1 s x s i with y (-1) = A and y(l) =B (1.22) 


where e is a small parameter (0 < c « 1 ) . he used a simple 

variable mesh approximation 


2e ^1+1+ Piyi-rfPi+q^Yi] , , 4 [«iy i+1 - p^i-i +<p r <3^] 

+ a(x^) 


p^itPi^i) 


p^itPi+qi) 
+ b(x i )y 1 = ffx.^ 


(1.23) 


where p.= x. +1 -x. and q ± = x. - x._ r 

This is a first order approximation, if we assume the mesh ratio 
gi/Pi as constant. A difficult solution procedure is followed, 
which finally requires a few thousand mesh points in the interval 
[-1, 1] . The mesh must be properly chosen so that the solution of 
the difference equation approximates that of the differential 
equation. This is accomplished by iteratively adjusting the mesh 
spacing such that the mesh points are concentrated in the region 

where y(x) changes rapidly. The procedure is first to solve the 

. . -2 
problem with a uniform mesh and a modest value of e ( c = 10 or 

—3 . . . 

10 ) . Then insert new mesh points between adjacent points, say x^ 

and for which a certain predetermined tolerance is exceeded 

(e.g., | Yji - y^ | > 5, 5 prescribed). A smoothing process is 

carried out to avoid locally abrupt changes in the mesh intervals. 

Next the problem is solved by Gaussian elimination. Finally, the 

size of e is decreased and the procedure is repeated, using the 
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mesh obtained from the preceding e-step as the initial mesh. 
Obviously, this method is costly in terms of computer time even for 
simple linear problems. 

Pearson [1968b] has also given the solution of the nonlinear 
problem 

F ( x ,y#y / , e y' 7 ) = o for o * x =* 1 ( 1 . 24 ) 

with y (0) = a and y(l) = £ (1.25) 


The same approximation as in case of linear problems is used for y' 
and y' ' . The resulting non-linear algebraic equations were solved 
by the Newton-Raphson iterative scheme. For example, to solve the 
non-linear problem 

e y" + y' 2 = 1 (1.26) 

y(o) = y ( i) = i (i.27) 

—6 

where e = 10 , Pearson used 4000 mesh points to get a solution 

which agreed up to five significant digits with the exact solution. 
In a more difficult problem, he used 25000 points to get an 
accurate solution. In order to overcome the numerical instability 
of the standard methods the upstream (upwind) one-sided 
(directional) differences on a uniform mesh is introduced. For 
detailed discussion see Dorr(1970). The essential of this type of 
scheme is to replace the first order derivatives in (1.21) by a 
one-sided difference quotient instead of the central difference. 
The choice of a backward or forward difference depends on the sign 
of a(x) at the particular mesh point under consideration, that is 


'a(x i )(y i+1 - y^/h if a (x i ) a 0 
3(5^)^ -yi_i)/h if a (x i ) < 0 
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Then the difference scheme takes the form 

r y i+i" 2y i + y i-i 


VYil - 


+ a [ ay^] + = f ± (1.29) 


If b(x^) ^ 0 then it can been shown that is a difference 

operator of the positive type and hence there exits a unique 
solution for each set of given data and for each c > 0, h > 0. By a 
difference operator of positive type we mean a operator of the 
form 


VYi> - A i y i-1 + B i y i + C i y i + 1 

where 

(i) A_ l > 0, V i, and (ii) 0, V i. 

With these above restrictions on b(x), satisfies a discrete 

maximum principle. Moreover, if the directional difference method 
is used, then for fixed h > 0 there is a limit function, as e tends 
to zero, which satisfies the 'reduced' difference equation, (1.29) 
with e = 0. One drawback of this method is that it is only of first 
order . 

Dorr et al. (1973) have developed the idea further and they 
discussed the applications of the maximum principle to obtain 
elementary estimates for solution of second order ordinary 
differential equations. These estimates are applied to obtain 
results on the limiting behavior of solutions of singularly 
perturbed problems. They considered the linear problems under 
various hypothesis, including turning point problems, and 
qausilinear problems of the form: 
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ey' 7 + a(t,y(t),e)y 7 (t) = 0 (t,y (t) ,e) 
for a < t < b with y(a) = A(e) and y(b) = B(c) 


(1.31) 

(1.32) 


II 'in (1969) constructed a finite difference scheme which represent 
the rate of decay in the boundary layer for the equation of the 
form: 

ey 77 + a(x) y 7 = f(x), (1.33) 
where a(x) s a is constant, 

and then it is generalised for slightly more general problem i.e., 
when a(x) is a variable and for a uniform mesh, the following 
difference operator is obtained, 


= T ± [ 


1 + 1 


” 2y ± + y 


i+1 


] + a ( x i> [ 


^i+i- y i-i 
2h 


]- 


f(x.) 


(1.34) 

. 2 
where y^ is chosen so that the scheme is 0(h) accurate and 

correctly represent the rate of decay from the boundary layer into 

the interior. He found that for the eq(1.33), 


a(x.)h r a(x.)h 

= 2 ~ coth 1^ ^ J (1-35) 

and so the operator is of positive type. 

. 2 

He then proved the if a and f are C on the interval, then 

I y^) - | * K h 2 , (1.36) 

where K is a constant which depends on e . 

Abrahamsson et. al.(1974) have also given finite difference scheme 
for singular perturbation problems where the equation is a system 


of the form: 
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ey" + A(x) y' + B(x) y = F(x) 0 s x s 1 (1.37) 

with boundary conditions 

y(0) = a and y(l) = £ (1.38) 

where e > 0 is a small parameter, y = (y (1) ,y (2) , • • *y (n) ) T and F = 
(F {1) , f ( 2) , . . . F (n) ) T are vector functions with n components and A,B 
e C°° are n-rowed square matrices. They assumed that 


A (x) 


A I (x)0 
0 A 11 (x) 


* . . I 

with A = A (=ad 3 omt of A) , A s -7) < 0, 


II 


== 7) > 0 for 0 ^ x 


They devised some schemes which with net spacing h » e, yield 
'accurate' solutions in the interior, that is away from boundary 
layers. This is done and Richardson Extrapolation is justified 
under appropriate assumption on A(x) but in general the accuracy 
cannot be better than 0(e) . In the scalar case, their work can be 
considered as refinement of upstream one-sided difference scheme. 
The idea is to introduce a parameter into the difference equation 
in such a way that more accurate approximation for the reduced 
problem is obtained. 

Nijima(1978) derived a difference scheme for a semilinear problem 
with any number of turning points of arbitrary orders. The class of 
the problems considered are 

ey" - [a(x)y]' - b(x,y) =0; 0 ^ x ^ 1 (1.39) 

with y(0) = A and y(l)=B (1.40) 

They showed that a solution of their scheme converges uniformly in 
c, to that of continuous problem. 



Lorenz (1979) has considered the numerical solution of singularly 
perturbed two-point boundary value problem of the form : 
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-ey" + a(x,y)y' + 0(x,y) = 0 (1.41) 

for 0 3 x * 1 with y (0) = A and y(l) = B (1.42) 

He essentially breaks the interval into two parts and uses 
different meshes on each part to solve a reduced and boundary layer 
problem. It is noted that the criteria for choosing the break point 
is not given although the cutting parameter cannot be too large. 
Further the solution of the reduced problem is used as the data for 
the boundary layer problem at the cut point. Thus Lorenz solves the 
reduced problem only upto the break point and then solves the 
boundary layer problem for the rest of the interval. 

Kriess and Kriess(1981) have given some finite difference schemes 
for a system of the form : 

^ = A(x) y + F (x) 0 s x s l (1.43) 

with n linearly independent boundary conditions 

R 0 y(0) + R x y(l) = g (1.44) 

where y = (y (1) ,y (2) , . . .y (n) ) T is a vector function with n 

components and R Q , R 1 and A(x) are n - rowed square matrices. 
Osher(1981) has developed upwind finite difference approximation 
for non-linear problems and also for system of non-linear 
hyperbolic conservation laws. Similar discussions with one sided 
difference schemes for the problem 

ey" - a(y)' - b(x,y) = F(x); -1 * x * 1 


with 


y(-l) = A and y(l) = B 


(1.45) 

(1.46) 
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where A,B are constants, a(y),b(x,y) are C 2 functions with 

b(x,0) = 0 and ^ 0 is contained in 0sher(1981). 

Abrahamsson and Osher(1982) have discussed the monotone difference 
schemes for the problem 

ey" - [f(y)]' - b(x # y) = 0, 0 3 x S 1 (1.47) 

with y (0) = A and y(l) = B (1.48) 


where f e C 1 (IR) and b e C 1 (IR 2 ) and > s > 0. 

They proved some convergence results and showed that the 
Engquist-Osher (1981) monotone scheme will reproduce essential 
properties of the true solution for any grid. 

Hoppensteadt and Miranker (1983) have discussed an extrapolation 
method for the numerical solution of singular perturbation initial 
value problems of the form: 

gf = f(t,x,y,e); x(0) = £(e) 


g| = g(t,x,y,e) ; y(0) = tj(e) 


(1.49) 


Sakai ( 1975) and Sakai and Usmani (1984) have discussed some 

numerical method based on chopping procedures for the solution of 
singularly perturbed two-point boundary value problems. 
Niijima(1984) has derived a completely exponentially fitted 
difference scheme for turning point problems. 

O'Malley and Flaherty (1980) have discussed some analytical and 
numerical methods for nonlinear singular singularly perturbed 
initial value problems. These results are extension of Ascher and 
Weiss (1983 , 1984a, 1984b) . Singular-singularly perturbed boundary 
value problems are also considered by O'Malley (1979). 
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Roberts ( 1982 ) has given boundary value technique to solve singular 
perturbation problems. It is based on the shooting methods. 

Roberts ( 1983 ) has also discussed the analytic and approximate 
solutions of the problem: 

ey" = yy' , -l * x * l ( 1 . 50 ) 

y(-l) = a and y(l) = /3 (1.51) 

Roberts (1984) has extended his boundary value technique to solve 
the problem: 

ey" + yy' - y = o; o =s x s i (1.52) 

y(0) = a and y(l) = 0 (1.53) 

Jain et al. (1984) have given a family of third order variable mesh 
difference methods for two-point, second order, singular 
perturbation problems of the form: 

Y" = f(x,y,y') (1.54) 

where f contains the parameter implicitly, Under the natural, mixed 
and nonlinear boundary conditions, when the mesh ratio is equal to 
unity, then these third order variable mesh method reduces to a 
unique fourth order method. 

Vulanovic(1989) has considered finite difference scheme for the 
following singularly perturbed boundary value problem: 

Tu s -ey" + b(u)u' + c(x,u) =0, x e I = [0,1] 

Bu = (u(0),u(l)) = (U Q ,U 1 ) 

Under the assumptions; 

HI. Let k = 1 or K = 2 and let 

b e C^(W) , c e C^(I x W) , where W = [u^,u*]; 

_ * , 

and u^ < u satisfy: 


(1.55) 

(1.56) 
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H2 . C(x,u*) 2: 0 a C(x,uJ, x e I 
U* ^ U. 2 U*, j = 0,1 
and C u (x,u) 2: c* > 0, x € I,u e W 

Vulanovic(1989) has also given a numerical method for quasilinear 
singular perturbation problems of the following form: 

Tu s -eu" + b(x,u)u 7 + c (x,u) =0, x e I = [0,1] (1.57) 

Bu= (u(0) ,u(l) ) = (Uq,^) (1.58) 

under the assumptions, 

b,c e C 2 (IxIR) , 

b(x,u) 2 : p > o, b x (X,U) + C u (X,U) 2 r, r S o, X € I, U € IR, 

£ 2 + 4 £K > 0 

The continuous problem is transformed to the form: 

eu" + f(x,u) 7 = g(x,u) , (1.59) 

where 

f(x,u) = S b(x, s) ds, g(x,u) = f x (x,u) + c(x,u) 

Then new variables (t,y) , are introduced which are given by 
x = X(t), Y = u(X(t)), where X e C 2 (I), X 7 (t) > 0 for tel and 

X (0) = 0, X ( 1) = 1. As a result of transformation following 
equation is obtained: 

Py a -e(M(t)y / ) / - p(t,y) 7 + q(t,y) = 0, (1.60) 

By = (Uq,!^) (1.61) 

where n(t) = 1/X 7 (t) , p(t,y) = f (X(t) ,y) ,q(t,y) = g(X (t) ,y) X 7 (t) . 
Then a finite difference scheme is given and first order 
convergence in the perturbation parameter is proved in discrete 

r 1 

L -norm. 
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Herceg(1990) has developed a uniform fourth order difference scheme 
for the singularly perturbed problem of the form: 

T c ~ ~ £ 2 u 7 ' + c(x,u) =0, x € i = [0,1], (1.62) 

u(0) = u(l) = 0 (1.63) 

under the following assumptions: 

c e C k (£ x 1R) , k € N, 

g(x) c u (x,u) =s G(x), (x , u) € l X !R, 

5 = min { 5g(x) - 2G(x) : x e l } > 0, 
o < K 2 * g(x), I g' (x) I S L, I G'(x) | ^ L, X € £. 

1.4.2 FINITE ELEMENT METHODS : 

Zienkiewicz et al. (1975) were the first to treat singular 
perturbation problems by finite element method. They have stated 
that a finite element equivalent to upwind difference was needed to 
avoid the problems of oscillatory solutions obtained with standard 
numerical approximations. A further algorithm using piecewise 
linear triangular elements in two dimensions was proposed by 
Tabata(1977) and a further evaluation of the methods as well as a 
comparison with equivalent finite difference formulations for the 
convective-diffusion equation has been given in Zienkiewicz and 
Heinrich (1978) . Some application of theory of Babuska (1973 , 1981) to 
the singular perturbation problems have been published by 
Reinhardt (1982) . 

Major contribution to the application of finite element methods for 
singular perturbation problems have come from Mitchell, 
Zeinkiewicz, Babuska, Hemker, Miller, Griffiths, Van Veldhuizen, 
Reinhardt and Christie. 
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Miller (1975) has given the concept of introducing a parameter into 
the numerical scheme and then choosing it in order to meet some 
criterion for singular perturbation problems of the form: 

-eu" + a.jU' + a Q u = f in Q = (0,1) (1.64) 
with u(0) = u (1) = 0 (1.65) 


where a Q ^ 0 and a^ £ 0 are constants. He treats this problem by 
means of finite element method. The parameter 0 = 9(c) is chosen so 


that 


a 1 0h 

0 3 0 £ 1, limit 9 = 0, limit - ■ = 1 


e — > 0 

where h is the uniform mesh size. 


2e 


( 1 . 66 ) 


N-l 


This parameter 9 is introduced in the basis as follows: 

I" (x-Xj_ 1 )/eh for x e [Xj_ 1 ,Xj_ 1 + eh] 

1 for x e [Xj_ 1 +eh, Xj] 


<P j(x) = 


(Xj+eh-x)/eh for x e [Xj,x^+ eh ] 


otherwise 


It is then shown that the finite element scheme is equivalent to 


e A a ll ^ r 

2e , , 


c a l 


eh 2 J u j+l 

0h + a o h J 

U . + 

3 

eh 2 

U j-1 


= h [f(x j _ 1 + eh) + f(x.)]/2, j = 1(1) N-l (1.67) 

with u Q = u N = 0 (1.68) 

In the limit when e = 0 and thus 0=0, because of (1.66), (1.67) 
becomes the upwind finite difference scheme for the reduced problem 
of (1.64), 
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u . — u . . 

* j 3’ 1 + a 
a l h + 0 ~j 


,) + f (x.) 
u . = J - . =L_ 


with 


U 0 = 0 


j=l (1) N-l (1.69) 

(1.70) 

In the limit when 9 = 1, (1.67) is the usual central difference for 

(1.64) 

u . . „ - u . 


-e 


u i +1 - 2U -, +u i-i 


+ a l ~ ^~ ~2h 3-1 + a 0 u j = f < x j> 


A particularly interesting intermediate choice of 0 is 


(1.71) 



r a i h -i 


r "1 

0 = 

[ tanh 2e J 

/ 

! 

1 ajh/ 2e J 


(1.72) 


Now, this 0 satisfies (1.66) and (1.67) becomes 


a l h 


coth 


(£P 


i+i- 2u i +u i-i 


+ SL. 


U j+l“ U j-l 

2h 


a 0 U D 


f(x. + 0h) + f(x.) 

= - i- j=l ( 1) N-l (1.73) 

with u Q =u N =0 (1.74) 

which is the finite difference scheme introduced by II 'in (1969) 
Miller has not carried out an error analysis in this paper, but has 
extended the above concept for partial differential equations in 
Miller ( 1976) . 

Christie et al. (1976) has considered the following problem: 

ey 77 - Ky 7 = 0, 0 < x < 1 (1-75) 

with u (0) = 1 and u(l) = 0 (1-76) 


He has placed a parameter in the test functions instead of using 
the parameter in the elements for obtaining asymmetric linear and 
quadratic basis functions for the test space. In this way the 
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oscillations which occur with symmetric test functions are no 
longer present. 

Veldhuizen(1978) has used a finite element method using piecewise 
polynomials of degree s k to approximate the problem 

eu" + u # « f; 0 s x =£ 1 (1.77) 

with u(0) = u(l) = 0 (1.78) 

with a very irregular mesh, on this mesh, error estimates of order 
kti 

0(h ) are obtained uniformly in e , h the maximum step size, for k 

^ 2 . They have extended these results and discussed some higher 
order schemes of positive type for the following class of problems 

eu" + a(x)u' = f(x) ; 0 ^ x ^ 1 (1,79) 

with u(0) = u(l) = 0 (1.80) 

where a(x) £ a Q > 0. 

They have also extended these results for elliptic singularly 
perturbed problems in two dimensions. 

Hsiao and Jordan (1979) have discussed numerical schemes based on 
the method of matched asymptotic expansion and modifying the 
boundary layer problem. The particular method they used for solving 
the modified problems in Galerkin method with linear finite 
elements as trial functions. They have applied these schemes 
successfully to several examples. Finally, an error analysis has 
been performed and encompasses both numerical and singular 
perturbation theory. 
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0 / Riordan(1984) has introduced a new piecewise linear finite 
element, which is designed to handle singular perturbation 
problems. They first introduced the new concept of hinged elements. 
Then they examined finite element approximations, to the problems 
of the form 

ey" + ay' = f; 0 s x s 1 (1.81) 

with y (0) = y (1) = 0 (1.82) 

whose errors are uniformly in e i.e. the error constants and 
asymptotic order of convergence as h — > 0 are independent of e . 
Piecewise exponential elements yield 

max | y (x. ) - y h (x. ) | s Ch 2 (1.83) 

0<i=sN 1 1 

II Y“Y h IL - Ch (1-84) 

where C is a constant independent of c and h. 

1.4.3 SPLINE APPROXIMATION METHODS : 

Starting from late sixties, many researchers have used splines 
for solving regular perturbation problems and some have also used 
it for solving singular perturbation problems. In this section, we 
give a brief survey of these methods. 

Loscalzo and Talbot(1967) were perhaps first to use spline for 
solving differential equations. They have given a procedure for 
obtaining spline function approximations for solutions of the 
initial value problems. He has used quadratic and cubic splines and 
showed that the method is related to well known trapezoidal rule 
and Milne-Simpson method and the method is shown to be divergent, 
when higher degree spline functions are used. 
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Bickley (1968) has used cubic splines for solving two-point boundary 
value problems. He has considered the general linear two point 
boundary value problems of the form: 


P(x)u" + q(x) u' + r (x) = s(x) 


(1.85) 


with boundary conditions 

oc Q u + p Q \i' = 7 q at x = x q (1.86) 

a n u " &n u ' = r n at x = x n (1 * 87) 
The essential feature of his analysis is that it leads to a set of 
linear equations whose matrix of coefficients is of upper 
Hessenberg form i.e. an upper triangular matrix with a single 
subdiagonal. 

Soon afterward, Albasing and Hoskin(1969) has used Ahlberg's cubic 
spline for obtaining the solution of the second order linear 
differential equation 


y" + f (x) y' + g(x)y = r(x) (1.88) 
with y(x Q ) = a and y ( x N^ “ b (1.89) 
He has considered separately the cases where first derivative term 
is absent and where the first derivative is present. 

For the first case, using the continuity conditions of first 
derivatives at the nodal points, he has obtained the following 
tridiagonal scheme, 


•)«( 1 * I Vi ) 




■5 2h rr 

2 - -3 g 


h 4 


i ) 1 * 5 


■ t ! ( 


j+1 + 4r j + r j“l J 


(1.90) 


and for the second case, he obtained the following scheme 
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( 1 “ 2 f j-*-i + « 9-i 


A j 


y j+H - 2 -j+1 ■ 6 3 j+l J 

- ^[f 1 + I f j+ i) A j +1 1 - ! f j+ i) B j - ^Vj - yj-ifi + |% j+1 ) 
j-i( 1 - I f j-i + iVi ) B j = | 2 ( A j r i+i + 4C j r j + Vj-1 ' 


+ y 


j = 1, 2, ... , n-1 


(1.91) 


where = 




1 - — f . 

3 r j-l J 


' 1 + ! f j ) + £ f 3-l f j 


B j A j+1 


C . = 1 + — 
j 24 


b 


j+i ■ f j-i j 


12 j-1 j+1 


Fyfe(1969) has examined the Bickley's method and obtained the error 
estimates. He also examined deferred corrections, effect of 
nonuniform spacing and a mesh refinement process. 

Jain, P. C. and Holla (1978) has given a scheme for linear 
hyperbolic equations with variable coefficients, where they have 
approximated the time derivatives by central differences and the 
space derivatives by second order derivatives of interpolating 
cubic spline. 

Patrico (1978) has presented a numerical process which provides 

4 . . . . 

0(h ) cubic spline approximate function approximation for the 

solution of general initial value problem. 

Patrico (1979) has presented a numerical process which provides 
0(h ) spline approximate function approximation for the solution of 
general initial value problem. He has used spline functions of 
degree five for obtaining the approximate solution. He has shown 
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that the method is stable and convergence analysis is given. 

Mastro and Voss (1981) have considered a quintic spline collocation 
procedure for solving the Falkner-Skan boundary layer equation. 

Jain and Aziz (1981) have constructed the following parametric spline 
2 

S " (t) = " w\in w r (t > )Sin W ( (t ~Vl )/h ) +S '' ( Vl )sin w (( t ? t)/h )] 

2 2 
+ [(( t -Vi)/ h )[s"(t l ,) +^S( V ) 

+ ((t-t)/h) (s"(V l)+ ^ S( Vl ))], 

(1.92) 

where w = h\/p, p > 0 and which reduces to cubic spline for p = 0. 

They have used this spline for solving certain ordinary 
differential equations and partial differential equations. They 
have compared results using parametric and cubic splines. 

Chawla and Subramanian(1987) have given a fourth order method based 
on Numerov's method and cubic spline for the following problem: 

Y" + f(x,y) =0, 0 s x * 1, (1.93) 

y (o) = <x, y (i) = £ (i.94) 

where f(x,y) is continuous, Sf/dy exists and is continuous 
First, They have extended Bickley's idea to the above non-linear 
problem and carried out analysis which shows that it is only second 
order accurate. Then they have given a fourth order method where 
Numerov's method is used to compute solution at the nodal points. 
These values of the solution are then used m continuity condition 
to obtain approximation of second order derivatives needed for 
construction of fourth order approximate cubic spline solutions. 
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They have presented numerical examples which demonstrate the fourth 
order accuracy of their method. 

Chawla(1988) has extended the results of above paper for general 
two-point boundary value problem of the form: 

y" = f(x,y,y'); a < x < b (1.95) 

subject to boundary conditions 

« 0 y(a) - o^y'fa) = A, /3 Q y(b) + fi ± y' (b) = B, (1.96) 

ot C)/ a 1 ,^ 1 and @ 2 are suitable constants. 

Surla and Stajonovic(1981) have give a difference 

scheme for singularly perturbed boundary value of the form: 

-e y" + p(x)y = f(x), p(x) > 0, (1.97) 

Y(0) = « 0 / Y(l) = (1.98) 

They have used the following spline in tension. 


>-j (x) = U.t + (l-t)Uj -;L 


+ Mj.! (h 2 /p 2 ) ( - (l-t)j, (1.99) 

t = (x-Xj^J/h, x e [ x j-i' x j3< Xj = jh, j = 0(1) n, h = 1/n, 

p = p h, p = (p/e) 1/2 

They have obtained error of the form 0(h,min(h,h/e) for the 
solution at the grid points and second order convergence for the 
solution at off-nodal points for the case p(x) = p > 0. 

Sakai and Usmani(1989) have considered an application of simple 
exponential spline to the numerical solution of the following 
singular perturbation problem: 
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ey" + b (x) y' - d(x)y = f(x) 0 =s x ^ 1 (1.100) 

Y (0) = a QI y (1) = /3 (1.101) 

for small e > 0 and for smooth data function b,d and f subject to 
the conditions d(x) 2 : 0 and b(x) i B > 0 on [0,1]. 

Stajonovic(1990) has considered the self-adjoint singular 
perturbation problem 


Lu 3 - cu"(x) + p (x) u (x) = f (x) , p (x) — p > 0 , x € (0,1) (1.102) 

u(0) = a Q , u(l) = a ±t cc 0 ,oc 1 € !R (1.103) 


and has used the following exponential spline difference scheme for 
obtaining numerical solution of the above problem: 


r I v i-i + r i v i + r i v i+ i - <*?!-! + 4 f i + 4 f i + i 


i=l(l) n-1 


v 0 a 0' 


V N " *1 


(1.104) 

(1.105) 


where , 

rT = A i _ 1 /2sinh(\ i _ 1 /2) , r+ = x i+1 /2sinh(A i+1 /2) , 
r? = - A.^coshX^/sinh(\^/2) 

q i = p i-l (1 “ A i-l /2sinh(A i-l /2)) / “ A i+i/ 2sinh (* i+1 /2) ) , 

q? = ( -2 + \ ^ cosh A./smh(A i /2)) 

P i = p(x i>' X i = (P i /e) 1 / 2 h. 

He has proved that for p(x), f(x) e C 2 (0,1] and p' (0) = p' (1) = 0 

| u(x^) - v^ | * M min(h 2 ,c) (1.106) 

where v^ denote the exponential spline solution and M is a 
constant . 


Roos(1990) has considered a self-ad 3 oint model singular 


perturbation problem 
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L^u s - e 2 u 7 ' + p(x)u = f(x), (1.107) 

u(0) = u(l) =0. (1.108) 

where it is assumed that p(x) £ p Q > 0 and p,f are continuous. 

He has constructed global uniformly first and second order schemes 
using patched base spline functions by replacing the variable 
coefficients of the differential equation by piecewise polynomials 
and to solve the resulting problem exactly. 

Bhatta and Sastri (1991) have given a seventh order spline 
procedure, where an octic spline coupled with a heptic spline 
function is used for solving the two-point boundary value of the 
form: 

y 7 7 = Py + Q (1.109) 

y (0) = a, y (1) = £ (1.110) 

where P and Q are continuous with P > 0, x € [0,1]. 


1.5 DEFINITIONS AND ALGORITHMS : 

1.5.1 Definition : A tria^uJra-r matrix M = ( m ij) of order n is said 

to be irreducible if 

m. . - * 0 2 ^ i - n 

l , ± 

and m^ ^ +1 * 0 1 * i ^ n-1 


1.5.2 Definition : A matrix M = (m^j) of order n is said to be 


diagonally dominant if 


n 


l * l s n 


( 1 . 111 ) 


£ i «ij i s i m i.i i' 
j=i 
j*i 

and strictly diagonally dominant if strict inequality holds 
for all i 
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1.5.3 Definition : A matrix M = ( m tj) order n is said 

ke irreducibly diagonally dominant if M is irreducible, 
diagonal dominant and strict inequality holds in (1.111) 
at least for one i. 

1.5.4 Definition : A matrix M = (m^j) of order n is said to be 
monotone if a o. 

1.5.5 Definition : If a matrix M = (m^) is irreducibly diagonally 
dominant and has nonpositive off-diagonal elements, then M is 
monotone. 


For solving the tridiagonal systems in this thesis, we have used 
the Thomas Algorithm which is described below : 


1.5.6 THOMAS ALGORITHM : Consider the following tridiagonal system 

of order n 


d l X l + °1 X 2 

a 2 x l + d 2 X 2 + C 2 X 3 

a_x 0 + d_x_ + 

3 2 3 3 3 4 



a _ i + d„ ,x„ , + c ,x n 
n-1 n-2 n-1 n-1 n-1 n 


= b 


n-1 


ax , + c x 
n n-1 n n 


= b 


n 


Assuming d 1 * 0, we eliminate x 1 from the second equation, getting 
the equation 


d 2 X 2 + C 2 X 3 = b 2 


= d 2 " d7 C l' 


b 2 - b 2 “ df b l' 


with d' 2 



Next, assuming * 0, we use this equation to eliminate x 2 from 
the third equation, getting the new equation 
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d 3 X 3 + C 3 X 4 = b 3 

. a 3 a 3 

with d 3 = d 3 “ 37 - c 2 / b 3 ^ b 3 ” d^ b 2 ' 

2 2 

Continuing in this manner, we eliminate at the step k, x^ from the 
equation k+1 (assuming that d/ * 0), getting the new equation 


d k+l x k+l + c k+l x k+2 b k+l 


with 


d k+l d k+l 


k+1 

d/_ c k 


b 7 . . = b. 


k+1 


b/_ 


'k+1 k+l d£. "k 


for k = 1,2,... ,n-l 


During back substitution, we first get, assuming d^ ^ 0, 

b 7 

n 

x = TT 
n d 7 
n 


and then, for k = n-l,...,2,l 


x k = 


b k - °k x k+i 
d k 


1.6 SUMMARY OF THE THESIS: 


In chapter II, we first, present a numerical method namely Method 
1, for Linear Singularly Perturbed Initial Value Problems. The 
method consists of dividing the inner and outer region and solving 
the inner region problem by a fourth order cubic spline after 
rescaling the inner region. The solution of the reduced problem 
provides the solution of the outer region problem. -Since The point 


of division of inner and outer region is arbitrary and also it is 
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observed that the solution obtained is discontinuous at the point 
of division. In an attempt to overcome these drawbacks, a more 
general method namely Method 2 is derived, which is applicable to 
Non-Linear Problems also. The method consists of a numerical 
cutting point technique. This method also takes care of continuity 
at the cut point. Both these methods provide global fourth order 
approximation to the exact solution and does not depend on 
asymptotic expansion. Computations are done and results are 
compared with exact solutions. An order table is presented which 
confirms the fourth order of the method. 

In chapter III, A class of Non-Linear Two-Point Singularly 
Perturbed Boundary Value Problems is Considered. The original 
Boundary value Problem is reduced to an asymptotically equivalent 
initial value problem by an initial value technique and is then 
solved in the inner region by cubic spline method for inner region 
employed in Method 2 developed in chapter 2 . The solution of the 
reduced problem of the original boundary value problem provides the 
solution in the outer region. Computations are done and the results 
are compared with solution obtained by others. Order table is 
presented for each example which gives the computational rate of 
convergence obtained by double mesh principle. 

In chapter IV, A class of general linear singularly perturbed 
boundary value problems is considered. We first consider the 
problem without first derivative term. Ahlberg's cubic spline in 
terms of moments is considered. A second order variable mesh finite 
difference scheme is developed using the continuity conditions of 
second derivatives of the spline at the interior nodes. Then the 
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variable mesh scheme is generalised for the problems with first 
derivative terms by taking some second order approximations. 

In Chapter V, A class of general semi-linear singularly perturbed 
boundary value problems is considered. We first considered the 
corresponding linear singularly perturbed problem. A family of 
variable mesh methods based on cubic spline approximations is 
developed, which gives third order approximations to the solution 
of the linear problem. Then these methods are used for solving 
semi-linear problems using a quasi-linearisation technique. When a 
parameter namely mesh ratio parameter present in the methods is 
taken to be unity, the family of third order methods reduces to a 
family of fourth order methods, and in addition, if other parameter 
introduced in the method is also taken to be unity, the family of 
methods reduces to a well known Numerov's method for the 
corresponding regular problem. 

In Chapter VI, the variable mesh methods based on cubic spline 
approximations developed in the previous chapter are generalised 
for the general nonlinear singularly perturbed boundary value 
problems . 

In each chapter, error analysis is given and computations are done 
by taking some well known examples. All computations are done on 
HP9000 supermini computer at IIT Kanpur in double precision. 



CHAPTER II 


Cubic Spline Methods for Singularly Perturbed 
Initial Value Problems 


2.1 INTRODUCTION : 

In many areas of application, notably fluid mechanics, 
electrical networks, chemical reactions and control theory, 
singularly perturbed initial value problems arise which have a 
narrow region known as 'initial layer' near the initial point 
where solution changes very rapidly. 

The use of conventional numerical methods for singularly 
perturbed initial value problems reguire a very fine mesh in the 
initial layer region, which makes these methods quite demanding on 
the computer time. Also the problem may become ill-posed 
numerically when mesh size gets too small. The use of asymptotic 
methods such as 'Matched Asymptotic Expansion' , which consists of 
(a) dividing the problem into an inner region problem and outer 
region problem, (b) expressing the inner and outer solutions as 
asymptotic expansions, (c) equating various terms in the inner and 
outer solutions to determine the various coefficients and (d) 
combining the inner and outer solutions in some fashion to obtain 
a uniformly valid solution, is not a routine exercise but requires 

skill and experimentation efforts. 

The object of this chapter is to present efficient numerical 
methods using cubic spline which give fourth order approximate 
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solution globally for singularly perturbed initial value problems 
(SPIVP) along with the smooth second order derivatives of the 
solution. These methods do not depend on asymptotic expansion and 
also do not require very fine mesh. 

First consider the Linear Singularly Perturbed Initial Value 
Problem (LSPIVP) 


ey' + p(x,e)y = q(x,e) a ^ x ^ b (2.1) 

y (a) = cx (2.2) 


where e is a small positive parameter, 0 < e « 1, « is a given 
constant, p(x,e) > 0 and q(x,e) are sufficiently smooth functions. 
Under these assumptions, (2.1) -(2.2) has a unique solution (Keller 
1968). Also it is expected that the solution of (2.1)-(2.2) would 
converge to the limiting solution u(x) as c » 0 where 


u(x) 


q(X/0 ) 

P(x, 0) 


(2.3) 


is the solution of the reduced problem at least away from an 
initial boundary layer region of non-uniform convergence. 

Then the method is modified and generalized for the 
Non-Linear Singularly Perturbed Initial Value Problem (NLSPIVP) . 

cy'= f(x,y,e) a * x * b (2.4) 

y(a) = « (2.5) 


Where is non-zero on (a,b) and has only one stable 

ay 

root for all values of x and y and there is an initial layer near 
x = a . Then the reduced problem 


f(x,u,0) = 0, 


( 2 . 6 ) 
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will provide the limiting solution u Q (x) and also it is expected 
that the solution of (2.4)-(2.5) would converge to the this 
solution of (2.6) as c — > 0 at least away from an initial boundary 
layer region of non-uniform convergence. Since Ug(a) * y(a) , the 
solution generally converges non-uniformly at x = a. 

In this chapter, numerical methods using a fourth order cubic 
spline for SPIVP's are presented. Our methods consists of: 

(a) Obtaining a solution of the reduced problem which provides the 
limiting solution; 

(b) Dividing the interval into an inner and an outer region; 

. . 4 

(c) Obtaining an 0(h ) approximation in the inner region using a 
cubic spline; 

(d) Obtaining a global approximation on the interval [a,b] by 

combining the solutions of inner and outer region problems. 

. . . 4 

These methods are analyzed and it is shown that they give 0(h ) 

convergence in the stretched variable over the inner region. 

2.2 METHOD 1 (For LSPIVP's) 

Using the stretching of the independent variable x by taking 
t = x/e, we transform LSPIVP (2.1) -(2. 2) to the following 
initial value problem 

gf + p(t,c ) 2 = q(t,e) (2.7) 

z(t Q ) = <*, (2.8) 

where t Q = a/e. 

2.2.1 CUBIC SPLINE METHOD FOR INNER REGION PROBLEM : 

Consider (2.7) -(2.8) on inner region [a/e, a/e + T] , where T 
is a positive integer. Given the values z Q ,z 1 ,...,z N of the 
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function z(t) at the equally spaced nodal points a/e = t Q , 
^1' * * * ' ^N-l' = a / e + T ' there exists an interpolating cubic 

spline S(t) with the following properties: 

(i) S(t) is a polynomial of degree three on [t. f t.,-];i = 0(1)N-1 

1 l+i 

(ii) S (t) is c 2 [t 0 ,t N ]. 

(iii) S (t i ) = z i# i= 0 ( 1) N. 

Consider the cubic spline given by (Ahlberg et al. 1967) 


S(t) = nu 


<*1+1-** <*-*!> <*-*!>■ (*1+1-*) 


- m. , , 
l+l 


(t i+l -t) 2 C 2 (t-ti) +h] (t-t i ) 2 [2 (t i+1 -t) +h] 

+ 2 . + * i+1 - ^3 


where nu - S' (t^) 


t L s t s t ±+1 , i « 0(1) N-l 


(2.9) 


and the second derivative is given by 


S" (t) = -2m i 


2t i+l + fc i ~ 3t 


- 2m i+l 


t. , , + 2t. - 3t 
1 + 1 1 


Z i+1 - Z i 


<*1+1 + t. - 2t ) . 


( 2 . 10 ) 


Now, the limiting values from the two sides are 


2m i 4m. z. - z. 

S"(t .-) = — + — - 6 , 

1 h h h z 


( 2 . 11 ) 


S"(t.+) = 


4m i 2m i+l , , z i+l' z i 
+ 6 g 

h h h z 


( 2 . 12 ) 
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The continuity condition imposed on interior nodes 
t^(i=l, . . . ,N-1) gives 

m i-l + 4m i + m i+i = E (z i+i " z i-i> ' i= 1 < 1 ) N_1 - < 2 - 13 ) 

Putting 

“i - Zj_ “ -P(t i( c) Zi + q(t i( E) i = 0 ( 1) N (2.14) 

in Eq.(2.13), we get the following two-step scheme for problem 
(2-7) — (2-8) : 

<E " p i-i )z i-i " 4 Pi z i“ <E + p i+i )z i+i " -^i-i + 4< 3i + « x+1 > 

1 = 1(1) N—1 (2.15) 

with z Q = a. 

which is fourth order Milne's scheme for the linear initial value 
problem. 

We next present our cubic spline method for solving (2.7)-(2.8) on 

~ 11 ~ in- 
step 1 : Using the above scheme, approximate solutions 

•fc * 

z 0 # z l' * * * /Z N f noc ^ a ^ Points are obtained. 

Step 2 : With z o' z i' * * * ' Z N' obtained in step 1, we compute 

* * 

a ® 

m* = -P(t i( e)z* + q(t i( e) i=l(l)N (2.16) 

i »0 = -P(t 0 ,c)a + q(t 0 ,c) 


with 
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Step 3 : With the values of z*,z*,...,z* and m* ,m*, . . . ,m* as 

0' 1' ' N 01 N 

obtained above, we construct the cubic spline S*(t) for 
(2.7) — (2.8) on t Q 3 t s t N by 


* * 


* * 
S (t) = 


(W* * (t-^) 2 (t i+1 -t) 


m i+l 


* (t i+r t} t 2 (t-t.)+h] 

+ z i rs + z i+i 


(t-t i ) 2 [2(t. +1 -t)+h] 


fc i “ " fc i+l' 1 = ° ( !) N-l 


(2.17) 


2.2.2 REDUCED PROBLEM : 


Solve (2.3) to get the solution u Q (x) which is the solution of 

(2.1) -(2. 2) in the outer region. 

Now finally we have global approximate solution y (x) of 

(2.1) — (2.2) on [a,b] as 


* 

y (x) 


S*(x/e) 

u Q (x) 


for x e [a,Tc] 
for x € (Te,b] 


(2.18) 


Note: In this method, there are two things which need some 
modifications. First, the cutting point, i.e. the point of 
division of inner and outer region, is arbitrary and secondly, the 
solution is discontinuous at this point. In the next section we 
have taken care of these points and also generalized the method for 
NLSPIVP's. 
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2.3 METHOD 2 (For NLSPIVP's) 

Using the stretching of the independent variable x by 
talcing t = x/e, we transform the NLSPIVP (2.4) -(2.5) to the 
following initial value problem 

d 7 

= f (t, z (t) ,e) (2.19) 

z(t Q ) = « (2.20) 

where t Q -a / e . 

2.3.1 REDUCED PROBLEM : 

Solve the reduced problem (2.6) and select the stable root 
u Q (x) i.e., for which is strictly negative. Then this 

root will provide the limiting solution of (2 . 4) - (2 . 5) . 

2.3.2 CUTTING POINT TECHNIQUE : 

Let t =t„+l be the initial cutting point. 

p 0 

For a positive integer N (not large), let h = (t p -t 0 )/N, 
tj^ — tg+kh, k=0 (1) N , tjj — tp . 

Step 1: For the initial value problem (2.19)-(2.20) over the 

interval [t Q ,t ], obtain the approximate solution z^, i = 0(1)N 
u p 

at nodal points t k , k = 0(1)N, by classical fourth order Runge- 
Kutta method (In fact, any fourth order method can be applied). 

Step 2: If |z N *- u Q (t p )| ^ h 3 , (2.21) 

then we stop and take t to be the desired cutting point, else we 

increase t by 1 and N by 2N then repeat step 1 until (2.21) is 
P 

satisfied. 

Thus, we have now obtained the desired value of the cutting 
point t and an approximation to the solution at the nodal points 
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over the inner region [t Q ,t ]. Since this solution z?, i = 0(1)N, 

r *■ 

is only a crude approximation, to obtain the solution (global) at 
the non-nodal points along with the smooth second derivatives in 
the inner region, we do the following: 

2.3.3 CUBIC SPLINE METHOD FOR INNER REGION PROBLEM : 

We next present our cubic spline method for solving 

(2.19)-(2.20) on t s t s t . 

P 

it it it 

Step 1 : With z Q , z i # * * • • z n- 1' obtained in the section 2.3.2, we 

it it 

compute m^, . . . ,m N _ 1 from the 'conditions of continuity' 


m i-l 

+ 4m^ + 

m i+i = 

•J 

h 

- ( z i+i - z i-i 

* 

* 

* 

3 

•k 

£ 

1 

to 

+ 

i 

M 

II 

+ 

h 

<W ‘ Z N- 


i = 1(1) N-2 ' 

► 


( 2 . 22 ) 


j 


with 

m* = f(t 0 ,«,€), = f(t N ,u Q (t N ),E) 

Since, it is cumbersome to write each time a separate equation for 

it 

replacing z N by u Q (t N ), therefore, now onwards, in this section 
we take z* = u Q (t N ). 

it it it it it it 

Step 2 : With the values of z Q ,z 1 ,...,z N and m 0 ,m 1 ,...,m N as 

it 

obtained above, we construct the cubic spline S (t) for 
(2. 19) -(2.20) on t Q s t s t using (2.17) of Method 1. 

it 

Now finally we have global approximate solution y (x) of 
(2 . 4) - (2 . 5) on [a,b] as 
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& 

S (x/e) for x e [a,x p ] 

y(x) = ■ (2.23) 

u 0 (x) for x 6 [x p( b] 

where x = t e 
P P 

Note : Since S (x^/e) = u 0 (tp), the solution is continuous at the 
cutting point. In fact, its first derivative is also continuous at 
this point. 

2.4 ERROR ANALYSIS ; 

Let e(t) denote the error in the spline solution S*(t) for the 
solution of (2 . 7) — (2 . 8) and (2.19)-(2.20) (Since the error 
analysis is more or less common for both methods, it is given 
together) . 

i.e., e(t) = z(t)-S*(t) = [z (t) - S(t)] + [S(t) - S*(t)] (2.24) 

= e x (t) + e 2 (t) (2.25) 

where e^(t) is the error due to spline interpolation 

and e 2 (t) is the error due to descretization of the 

differential equation. 

In the following, we use for v e C[a,b], || v || = max |v(t)|. 

[a, b] 

It is well known that the cubic spline interpolation is fourth 

4 

order process (Hilderbrand 1956). Therefore, for z e C [t Q ,t p ], we 
have 

|| e 1 | w s C^h 4 , C x a constant (2.2 6) 
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Next, for e 2 (t), from (2.9) and (2.17) we get 


e 2 (t) = (m i -m i ) 




+ (» i+1 - » i+1 ) 


C t-t ± )-( t i+1 -t) 


+ — - 1 , 3 ^ 


+ (Z i+l- Z i+l> 


* (t-t.)^[2(t. +1 -t)+h] 


for t i t ±+1 , i=0 ( 1) N-l 


(2.27) 


Let Z (z^,...,z^j) , Z— (z^, . . . , Zj|) 

M = (m^, ... /E 1 ^) / M = (m^, • • • 


In the following, we use for X = (x 1 , . . . ,x N ) , 


X II = max |x. 
11 00 . 1 
l^i<N 


from (2.27) we have 


where, 


e 2 (t) | s | M-M | f x (t) + | Z-Z | f 2 (t) 

s II M-m\ ^(t) + II Z-Z* I. f 2 (t) (2.28) 


( t - t) 2 (t-t ) (t - t ) 2 (t -t) 

f,(t) = 5 — + ^ 


f 2 (t) = 


(t i+1 -t) 2 [2(t-t i )+h] + (t-t i ) 2 [2(t i+1 -t)+h] 


and 
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Clearly, maxima of f 1 ( t) and f 2 (t) is attained at t = 


t +t. 

i i+l 


Hence, max f n (t) = —r— and max f (t) = 1 
t x 4 t 2 

Therefore, 


e 2 I. * T 


+ I z-z* || M 


(2.29) 


Bound for || M-M*|| m (For Method 1) 


From (2.14) and (2.16), we have 



m.- m . = 

-p(t i( e) ( z^- z i ) 


1 1 

= p(t if E) | z^ z*| 


1 M_M *I„ * 

p 1 


where P = max p(t^,e). 
i 

Bound for || M-M*|| w (For Method 2) 

Writing (2.13) in the matrix-vector form, we have 

A M = (3/h) B Z - C (2.31) 


Where A and B are tridiagonal matrices given by 
A i = i th row of A = Trid [141] 

B i = i th row of B = Trid [-101] i = 

and C = ( m Q ,0, . . . ,0,m N ) is an N-column vector, 

where m Q = f(t 0 ,z 0 ,e) and m N = f(t N ,z N ,e) 


1 , 


,N 
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Similarly, writing (2.22) in the matrix-vector form, we have 

A M*= (3/h) B Z*- C* (2.32) 

. * , * * 
where C = ( m Q , 0, . . . , 0,1^) is an N-column vector. 

Subtracting (2.32) from (2.31) we have 

A (M - M*) = (3/h) B (Z - Z*) + (C* - C) (2.33) 

From (2.33) 

II * (3A) HA" 1 ! || B |j |Z-Z*| M + || A** ^ | | C-C*| (2.34) 

Note that ||A 1 || ^ 1/2 ,also ||B|| = 2, 

® » s A 

From the definition of m^, and the mean value theorem, 

We have 

I ®H - ■S I S d I Z N “ W I 

* d (I Z N ~ Z N I + I 4 " W I 

^ d (|| Z~Z*|| co + h 3 ) (using (2.21) (2.35) 

where d = | ^(t N ,z N ) |, z N lying between z N and z N 

• * 

Now, from the definition of C and C , we have 

I c " c *l = I I 

which implies 

I C-C*| 3 d (|| Z-Z*||j- h 3 ) (2.36) 

Therefore, we obtain 

|M - M* ||„ * ( l + f)|| Z-Z*|| m + | h 3 


(2.37) 



Since we are using fourth order schemes in both methods for 
approximating z(t) at nodal points we have 
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1 2 - z * L * S h4 

where C 2 is a constant independent of h. 
Therefore, from (2 .29) , (2. 30) and (2.37) 

« e 2 I. £ C 3 h4 + C 4 h5 


where C 3 



(For Method 1) 
(For Method 2) 


where C. 

4 



(For Method 1) 
(For Method 2) 


Finally, combining (2.26) and (2.39) we have 

|| e l o S C h 4 , C = Cj_ + C 3 

Thus, we have proved the following: 

4 

Theorem 2.1: Let z(t) e C [t ,t ], then our cubic spline 

\J y 

4 . . * 

provide an order h convergent approximation S (t) for the 
z(t) of the initial value problems (2.7) -(2.8) and (2.19) 
that is, 

I • II* * c h 4 

where C is a constant independent of h. 


(2.38) 


(2.39) 


(2.40) 

methods 
solution 
-( 2 . 20 ) ; 
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2.5 NUMERICAL EXAMPLES : 

Example 2.1 . (Artificial test problem of Lapidus and Seinfeld 1971, 
Doolan et.al 1980) 

ey # = g(x) + eg' (x) - y 0 < x ^ 1 

y (0) = 10, g(x) = 10 - (10 + x)e“ X 

Here f(x,u,0) = 0 =» u = g(x) and f (x,g(x),0) = -1 

which is stable. Thus, u 0 ( x ) = <?(*) provides the limiting 

solution. 

The exact solution is given by 

y(x) = g(x) + 10 e“ x ^ c 

The computational results obtained by Method 1 and Method 2 are 
presented in Tables 2.1 and 2.2 respectively, for different values 
of h and e. 

Example 2 . 2 (Edsberg 1976, Doolan et.al 1980) 

2 

ey'=~2y 0 s x * 1 

y (0) - 1 

Here f(x,u,0) = 0 =* u = 0 and fy(x,0,0) = 0, but u = 0 is unique 
Therefore it provides the limiting solution u Q (x) = 0. 

The exact solution is given by 

y(x) « 1 / [1 + 2 x/c ] 

The computational results obtained by Method 2 are presented in 
Table 2.3 for the different values of h and e. 



Example 2 . 3 


1 x ^ 2 
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e Y' = *y(l + Y) 1^X52 

Y(l) = -2 

Here f(x,u,0) = 0 =» u = 0, -1 and fy(x,0,0) = x which is unstable 
while f y (x,-l,0) = -x which is stable. 

Thus, u Q (x) = -1 provides the limiting solution. 

The exact solution is given by 


y(x) = ex P((x -1)/ 2s + log 2 ) 

1 - exp((x 2 -l)/ 2c + log 2) 


The computational results obtained by Method 2 are presented in 
Table 2.4 for different values of h and c. 

2.6 DISCUSSION : 

4 

We have presented numerical method which gives 0(h ) 
approximation to the solution of singularly perturbed initial 
value problems. In the method l(For LSPIVP's) , choice of 
the cutting point is arbitrary and the solution at the cutting 
point is discontinuous. In method 2 (For NLSPIVP's) which is 
generalization and modification of method 1, a numerical cutting 
point technique, which is iterative in nature, is given for 
selecting the cutting point. Also the continuity of the solution 
at cutting point is guaranteed. 

Theoretically, problem (2 . 19) - (2 . 20) can be solved on 
interval [a/e,b/c] by any initial value method, practically c 
being very small, this interval is quite large. It is observed 
from the computation that by the cutting point technique, the 
effective size of the inner region in the stretched variable 
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remain almost same with the deceasing value of c . 

Also Milne's method can be used for obtaining solution at 

nodal points in method 2 and m's can be obtained subsequently as 

in method 1, but for solution to be continuous at the cutting 

point we have to change solution at this point from z* to u 0 (t ), 

N o p 

causing the violation of continuity conditions. In order to make 
solution smooth in the inner region along with matching it with 
outer solution at the cutting point, we have to recompute m's. 
Therefore solution at the nodal points serves only as fourth order 
approximations, which can be provided by any other fourth order 
method also. We preferred to use classical Runge-Kutta method for 
this purpose. 

Use of conventional methods for singularly perturbed initial 
value problems by transferring e to the right hand side generally 
require the step size to satisfy h « e even in the region away 
from the boundary layer region and thus restriction on the step 
size makes these methods quite demanding on the computer time. Our 
method does not involve any constraint on the step size h. 

Three test examples have been solved by this method. Example 
no.l is solved by method 1 and method 2 for comparison. We have 
computed solutions at the middle of the nodal points and we have 
tabulated the maximum and minimum error in both inner and outer 
regions along with the absolute dif ference (ds) between computed 
inner and outer computed solution at the cutting point for 
different values of h and e, it is observed that ds is non zero 
when we apply method 1 while it is zero for method 2, which 
confirms the continuity of the solution at the cutting point for 
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method 2 . Also the rate of convergence of these approximations 
which confirms the fourth order of our cubic spline methods in the 
inner region is given in table 2.5. From the tables it is observed 
that present methods approximates the exact solution well in both 
inner and outer regions and also the error does not grow with 
the decreasing values of e. In fact by our methods, which gives a 
global approximation of the solution, approximation at any point of 
the interval can be obtained without any further computation. 

We define the computational order of convergence for two successive 


values of h in the usual way i.e, for h and h/2 with respect to 


errors 


I h and I h ' 2 : 


Order = 


[ lo< ? ( I max) " lQ 9 ( I max) ] 


log 2 


Heading in the tables are as follows: 

h Mesh size 

I Maximum error in the inner region 

max 

I . Minimum error in the inner region 

min 

0 Maximum error in the outer region 

max 

0 . Minimum error in the outer region 

min 

ds Absolute difference between computed inner and outer 

solution at the cutting point 


I l T 
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Table 2.1(a) 




Example 2.1, e - 

10” 3 (Method 1) 



h 


T = 2 

T = 5 

T = 

10 


*max 

4.28810 E-5 

4.28810 E-5 

4.28810 

E-5 


I . 
mm 

0 

max 

2.41803 E-5 

6.01878 E-7 

6.01891 

E-7 

1/5 

1.35335 

6.73794 E-2 

4.53400 

E-4 


0 . 
min 

0.0 

0.0 

0.0 



ds 

1.35339 

6.74890 E-2 

1.38547 

E-4 


1/10 


1/20 


^max 

2.90740 

E-6 

2.90740 

E-6 

2.90740 

E-6 

^min 

1.86229 

E-6 

8.27880 

E-2 

1.18289 

E-10 

0 

max 

1.35335 


6.73794 

E-2 

4.53400 

E-4 

0 . 
min 

0.0 


0.0 


0.0 


ds 

1.35339 


6.73756 

E-2 

4.34932 

E-4 

■^max 

1.92246 

E-7 

1.92246 

E-7 

1.92246 

E-7 

nmin 

1.16134 

E-7 

1.12907 

E-8 

1.61059 

E-ll 

0 

max 

1.35335 


6.73794 

E-2 

4.53400 

E-4 

0 . 
min 

0.0 


o 

• 

o 


0.0 


ds 

1.35339 


6.73756 

E-2 

4.34932 

E-4 
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Table 2.1(b) 





Example 2 . : 

1, e =* 10” 5 (Method 1) 



h 


T = 

2 T = 5 

T = 

10 


■^max 

4.28810 

E-5 4.28810 E-5 

4.28810 

E-5 


^min 

2.41803 

E-5 6.01878 E-7 

6.01891 

E-7 

1/5 


0 

max 

0 . 
min 

1.35335 

0.0 

6.73794 E-2 

0.0 

4.53400 

0.0 

E-4 


ds 

1.35339 

6.74890 E-2 

1.38547 

E-4 


■^max 

2.90740 

E-6 2.90740 E-6 

2.90740 

E-6 


*min 

1.86229 

E-6 8.27880 E-2 

1.18289 

E-10 

1/10 


0 

max 

0 . 
mm 

1.35335 

0.0 

6.73794 E-2 

0.0 

4.53400 

0.0 

E-4 


ds 

1.35339 

6.73756 E-2 

4.34932 

E-4 


^max 

1.92246 

E-7 1.92246 E-7 

1.92246 

E-7 


^min 

1.16134 

E-7 1.12907 E-8 

1.61059 

E-ll 

1/20 _ . 


0 

max 

0 . 
min 

1.35335 

0.0 

6.73794 E-2 

0.0 

4.53400 

0.0 

E-4 


ds 

1.35339 

6.73756 E-2 

4.34932 

E-4 
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Table 2 

.2 



Example 2.1 

(Method 2) 

h 


e = 10~ 3 

x = 2E-2 

P 

-5 

e = 10 

X s ® 2E-4 

P 

1/5 

I max 

■^min 

0 

max 

0 . 
min 

ds 

4.45836 E -5 

6.32425 E -11 

2.06115 E -8 

0.0 

0.0 

4.45836 E -5 

7.82276 E -11 

2.06115 E -8 

0.0 

0.0 


J max 

2.49986 E -6 

2.49986 E -6 


^min 

3.41245 E -11 

4.27864 E -11 

1/10 


0 

max 

2.06115 E -8 

2.06115 E -8 


0 . 
min 

0.0 

o 

• 

o 


ds 

0.0 

0.0 


■^max 

1.52456 E -7 

1.52456 E -7 


^min 

1.45095 E -12 

2.49720 E -12 

1/20 


. 0 

max 

2.06115 E -8 

2.06115 E -8 


0 . 
min 

0.0 

0.0 


ds 

0.0 

0.0 
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Table 2 . 3 




Example 2 . 2 

(Method 2) 

h 


-3 

e = 10 

-5 

c = 10 



x = 1.2E-2 

P 

X = 1.2E-4 

P 


■^max 

8.65242 E -4 

8.65242 E -4 


■^min 

0 

max 

3.48361 E -9 

3.48361 E -9 

1/5 

2.73582 E -6 

2.73582 E -6 


0 . 
min 

o 

• 

o 

o 

• 

o 


ds 

0.0 

0.0 


"^max 

7.46802 E —5 

7.46802 E -5 

1/10 

"Snin 

0 

max 

3.50942 E -10 

2.73582 E -6 

3.50942 E -10 

2.73582 E -6 


0 

min 

o 

• 

o 

0.0 


ds 

0.0 

0.0 


"^max 

5.48070 E -6 

5.48070 E -6 

1/20 

■^min 

0 

max 

2.34987 E -11 

2.73582 E -6 

2.34987 E -11 

2.73582 E -6 


0 . 
mm 

o 

• 

o 

o 

• 

o 


ds 

o 

• 

o 

o 

• 

o 
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Table 2.4 




Example 2 . 3 

(Method 2) 


h 


-3 

e = 10 

e = 10~* 5 




x = 1+13E-2 

P 

x = 1+13E-4 

P 



^max 

3.89349 E -4 

3.89755 E -4 



^min 

0 

max 

0.0 

0.0 


1/5 

1.12921 E -6 

1.12921 E -6 



0 . 
min 

o 

• 

o 

o 

• 

o 



ds 

o 

• 

o 

o 

• 

o 



I max 

3.14718 E -4 

3.15041 E -4 


1/10 

"Snin 

0 

max 

0.0 

1.129211 E -6 

0.0 

1.12921 E -6 



0 . 
min 

0.0 

o 

• 

o 



ds 

o 

• 

o 

0.0 



^max 

2.21353 E -4 

2.21556 E -4 


1/20 

"^min 

• 

0 

max 

0.0 

1.129211 E -6 

0.0 

1.12921 E ~6 



0 . 
mm 

0.0 

0.0 



ds 

0.0 

0.0 
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Table 2.5 


h 

I _ 

Order 


max 

1/5 

4.28810 E -5 



Example 2 . 1 

1/10 

2.90740 

E -6 

3.88 

(Method 1) 

1/20 

1.92246 

E -7 

3.92 


1/5 

4.45836 

E -5 


Example 2 . 1 

1/10 

2.49985 

E -6 

4.16 

(Method 2) 

1/20 

1.52456 

E -7 

4.03 


1/5 

8.65242 

E -4 


Example 2 . 2 

1/10 

7.46802 

E -5 

3.53 

(Method 2) 

1/20 

5.48070 

E -6 

3.77 


1/5 

3.89349 

E -4 


Example 2 . 3 

1/10 

3.14718 

E -5 

3.63 

(Method 2) 

1/20 

2.21355 

E -6 

3.83 



CHAPTER ill 


On a Class of Non-Linear Singularly Perturbed 
Boundary Value Problems Via Initial Value Methods 

3. 1 INTRODUCTION : 

In general, finding numerical solution of a boundary value 
problem is more difficult than finding numerical solution of 
corresponding initial value problem. Therefore it is better to 
convert the second order problem into an asymptotically equivalent 
first order problem, wherever possible. Solving non-linear boundary 
value problems directly by finite difference schemes requires 
linearization. An iterative process with a fine initial guess and 
good number of iterations is then required to get an acceptable 
solution. 

In this chapter, we have given an alternative approach for 
solving a class of non-linear singularly perturbed boundary value 
problems, so that the above complexities can be avoided. Solving 
these problems by asymptotic methods is also far from trivial 
since it requires lot of experimental skill. We have converted a 
class of singularly perturbed non-linear boundary value problems 
into asymptotically equivalent singularly perturbed initial value 
problems. The cutting point technique is then applied to separate 
the inner and outer regions. The solution in the inner region is 
globally approximated by fourth order cubic spline method for 
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singularly perturbed initial value method, developed in the 
previous chapter . The solution of the reduced problem is taken as 
the solution in the outer region 

Consider a class of non-linear singular perturbation problem 
of the form 

ey , '(x) + [p(y (x) ) ] ' + q(x,y(x)) = r(x) a =s x * b (3.1) 

Y (a) = a, y (b) = j3; (3.2) 

here, e is a small parameter, 0 < e « 1; a and /3 are given 

constants; p(y), q(x,y) , r(x) are assumed to be sufficiently 

differentiable functions. Furthermore, we assume that the problem 
(3.1) -(3. 2) has a solution which displays a boundary layer of width 
o(e) at x = a for small value of e. 

This problem has been considered earlier by Kadalbajoo and 
Reddy (1988) . After reducing the original problem to an 
asymptotically equivalent singularly perturbed initial value 
problem, they have applied classical fourth order Runge-Kutta 
method. But since the resulting problem is still singularly 
perturbed, applying any conventional method for solving it, 
requires a very fine mesh (h « e) to satisfy the criterion of 

absolute stability for obtaining an acceptable approximate 

solution. In view of e being very small in practice, the 
restriction on the step size h would often lead to prohibitive cost 
of computational time. 
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3.2 INITIAL VALUE METHOD : 

We describe the technique in a step-wise manner as follows: 

Step 1. Obtain the reduced problem by setting e = 0 in (3.1) and 
solve it for the solution. Let y Q (x) be the solution of 
the reduced problem of (3-1) — (3.2) , that is 

[P(Y 0 (x) ) ]' + q(x,y o (x)) = r(x) (3.3) 

y o (b) = p (3.4) 

Step 2 Set up the approximate equation to Eq. (3.1) as follows: 

ey" (x) + [p(y(x) ) ]' + q(x,y Q (x)) = r(x) (3.5) 

where the y(x) term in q(x,y(x)) is replaced by y Q (x), the 
solution of the reduced problem (3.3)-(3.4). 

Step 3. Replace the approximated second-order problem (3.5) by an 
asymptotically equivalent first order problem as follows: 
By integrating (3.5), we obtain 

ey' (x) + p(y(x)) + Q (x) = R(x) + K (3.6) 

where 

Q(x) = fq(x,y 0 (x))dx, 

R(x) = jr(x)dx 

and K is a constant to be determined. 

In order to determine the constant K, the condition that the 
reduced equation of Eq.(3.6) should satisfy the boundary condition 
y(b) = /3 gives 
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P(y(b)) + Q(b) = R (b) + K (3.7) 

so that 

K - p(/3) + Q(b) - R(b) (3.8) 

Now, we adjoin the condition y(a) = a to (3.6) to obtain a 
first-order problem as follows: 

cy 7 (x) + p(y (x) ) + Q (x) = R(x) + K (3.9) 

y(a) = cl ( 3 . 10 ) 

where K is a constant given by Eq. (3.8). 

Step 4. Setting t = x/e in (3.9) and rescaling with 


y(x) = z(t) (3.11) 

y' (x) = z ' (t)/e (3.12) 

We obtain the resulting problem as 

z' (t) + p(z(t)) + Q (t) = R (t) + K (3.13) 

z(a/e) = a (3.14) 


Now, we apply the cutting point technique (Sec 2.3.2, Chap 2) for 
finding the cutting point t , then the cubic spline method for 
finding the solution of the inner region problem (Sec 2.3.3, Chap 2) 
is applied and solution y^ of the reduced problem is taken as the 
solution of the outer region problem. 

Step 5 Finally, the solutions of the inner and outer region 
we have global approximate solution 
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* S (x/e) for x e [a,x ] 

y W = \ p (3.15) 

l Y 0 (x) for x e [x »b] 

Jr 

where x^ = et^ and S (x/e) is the cubic spline solution obtained 
from Method 2 of second chapter. 

3.3 ERROR ANALYSIS ; 

Since, after reducing the original boundary value to an 
initial value problem and then applying the method 2 developed in 
the previous chapter, the error analysis of the present method is 
similar to that of above said method, therefore, in order to avoid 
repetition, we give only brief description of the error 
propagation. 

Let e(t) denote the error in the spline solution S*(t) for the 
solution of (3 . 13) - (3 . 14) , in the inner region [ a 0 / c ' t?] , 
i . e. , 

e(t) = z (t) -S* (t) - [ z ( t ) - S (t) ] + [S (t) - S* (t) ] (3.16) 

= e^t) + e 2 (t) (3.17) 

where e^t) is the error due to spline interpolation 

and e 2 (t) is the error due to descretization of the 

differential equation. 

Also, from Eq. (2.26) and (2.29) of Chap 2, we have 

|| e 1 || m ^ C^h 4 , C x a constant (3.18) 

II ® 2 I. a i 1 M " M * L + II z - z * IL < 3 - 19 > 
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where M = (n^, . . . ,m N ) t , M*= (m*, . . . ,m *) t , 

, .t _* . * * t 

Z - (z 1# ...,z N ) , z = (z 1# ...,z N ) are the actual and approximate 

solutions of (3.13) and (3.14) at the nodal points respectively. 
Also, 

C — (m Q , 0, . . . ,0,11^) , C = (m Q , 0, . . . ,0,1^) 


Where m Q = m* = -p(a) - Q(a Q /c) + R(a Q /e) + K, (3.20) 

m N = “P(z(t N )) - Q(t N ) + R(t N ) + K (3.21) 

m* = “P(Y 0 (t N )) - Q(t N ) + R(t N ) + K (3.22) 

Therefore, Using mean value theorem, we have 

I m N ” ”n I s d (1 Z ~ Z *IL + h3 ) < 3 - 23 ) 


where d = | §§(t N ,z H ) |, z N lying between z N and z* 

Now, from the definition of C and C , we have 

I C-C*| s d (|| Z-Z*^ + h 3 ) (3.24) 

Also from Eq(2.36) of Chap 2, we have 

1 M-m\ i (3/h) | A _1 1 1 Bj| 1 Z-Z* + IA _1 || | C-C* | (3.25) 

where, A and B are triangular matrices with ||A II - 2 and ll B ll = 2 
Therefore, we obtain 

» M - L s < E + f>» Z - Z *l» + I h3 (3 ' 26) 

Since we are using fourth order schemes in both methods for 
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approximating z(t) at nodal points we have 

a z - z* il * s h4 (3 - 27) 

where is a constant independent of h. 

Therefore, from (3.19), (3.26) and (3.27) 

» e 2 IL s C 3 h4 + C 4 h5 (3 - 28) 

7 ^ ^ 

where = 4 " ^2 + 8^ and — -g 

Finally, combining (3.18) and (3.28) we have 

a e IL s 5 h4 ' 5 = c i + c 2 (3 * 29) 

Thus, we have proved the following: 

Theorem 3.1 : Let z(t) e C 4 [a Q /e,t ], then our cubic spline method 
provides an order h 4 convergent approximation S (t) for the 
solution z(t) of the initial value problem (3 . 13) - (3 . 28) ; that is, 

a e il s 5 h4 

where C is a constant independent of h. 

3. 4 NUMERICAL EXAMPLES : 

Example 3.1 Consider the following example (O'Malley (1974), 
page 9, Eq. (1.10), Case 2) 

£y' * = yy ' — 1 ^ X — 1 

y (-1) = o, y(i) = -i 
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We have chosen O'Malley's approximate solution for comparison, 


Y(x) 


1 - exp f - 
1 + exp j - Mir 


For this example, we have a boundary layer of width 0(e) at 
x = -1. (O'Malley (1974) pages 9 & 10) 


Rewriting the differential equation in the form (3.1), we 

have 

ey" - ]'= o 


Integrating this, we get 
ey' (x) _ = K 


The constant K is determined using (3.8) as 


= . j 'hi) = 
2 


1 

2 


Hence the initial value problem is given by 


*»' - [ 
y (-i) - o 

This I VP is solved over the inner region m the stretched variable 
by the cubic spline method and the outer solution is given by the 
solution of the reduced problem 

Y q (x) Yq (x) = 0 
Y 0 (l) = -1 
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which gives Y Q (x) = -l. 

The computational results are presented in Tables 3.1 for different 
values of h and e. 


Example 3.2 Consider the following example (Kevorkian and Cole 
(1981) page 56, Eq. (2.5.1)): 

ey" + yy' - y = 0 O^x^l 

y (0) = -1, y (1) = 3.9995 

We have chosen to use Kevorkian and Cole's uniformly valid 
approximation for comparison: 


y(x) 


x + C 1 tanh 



where C^ = 2.9995, 

and C 2 = lo? [ C^+T 


For this example we have a boundary layer of width 0(e) at x = 0. 
(Kevorkian and cole 1981) page 56-66. 

Rewriting the differential equation in the form of (3.1), we have 


ey 


/ / 



= 0 


The reduced equation is given by 



This implies two equations 


y (x) = 0 
y' (x) = 1 
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Second equation is the correct differential equation, 
first equation is not satisfied at the boundary point x = 1 . 

Hence the reduced problem is given by 

Yq( x ) = 1 
y Q (x) = 3.9995 

whose solution is 

y Q (x) = x + 2.9995 
We get the approximate equation as 

ey" (x) + [ (x + 2.9995) = 0 

Now by integrating this, we get 

ey' + [ ( f~ + ^ -1 ) x ) “ K 

where £ = 3.9995 

The constant K is determined using (3.8) as 

K = ( | + (3 -1) = (g ~ 1) 

Hence the initial value problem is given by 
cy' (x) + [ Y(x) 2 - (x+g-1) 2 j _ Q 

y (o) = -i 


since 


This IVP is solved in the inner region in the stretched variable 
by the cubic spline method. 
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The computational results are presented in Table 3,2 for different 
values of h and e . 

Example 3.3 ; Consider the following example [Bender and Orszag 
(1978) page 463, Eq. (9.7.1)]: 

ey" + 2y' + e Y = 0 0 s x s l 

y (0) = 0, y(l) = 0 

We have chosen to use Bender and Orszag' s uniformly valid 
approximation for comparison. 

y(x) = log J- exp log2 

For this example we have a boundary layer at x = 0. 

The reduced problem is given by 

2Yq(x) + exp[y Q (x)] =0 0 s x s 1 

y 0 (D = 0 

which gives 

y 0 w = lo * 

We get the approximate equation as 

ey" (x) + [2y(x) ]' + = 0 

Now by integrating this, we have 

ey 7 (x) + 2y(x) + 2 log(l+x) = K 
The constant K is determined using (3.8) as 
K = 2y(l) +2 log (1+1) = 2 log2 
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Hence the initial value is given by 

cy' (x) + 2 |y(x) - log | ^ - ] 1 = o 

y (0) = o 

This IVP is solved in the inner region in the stretched variable 
by the cubic spline and the outer solution is given by the 
solution of the reduced problem. 

The computational results are presented in Tables 3.3 for different 
values of h and e. 

3 . 5 DISCUSSION : 

We have implemented the present method on three examples by 
taking different values of e and have compared the computational 
results with the approximate solutions developed by others. It can 
be observed that our method does not require a fine mesh and also 
does not impose any restriction on the step size h. In fact, 
numerical solution obtained for a crude mesh size h = 1/5 compare 
very well with the approximate solutions. Our method provides an 
alternative approach to the conventional way of solving certain 
class of non-linear singularly perturbed problems. The method does 
not require any asymptotic analysis and does not require any non 
linear equations to be solved. 

Since the exact solutions of the test examples are not available, 
we used the double mesh principle defined below for finding the 
computational rate of convergence: 
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Let 


z h = max I - u^ 2 | , 


j = 0 , . . . , N-l 


h . 


where u^ is the computed solution at the point t = t^ + Ah, 

A e [0,1], on the mesh ( t. \ , where t. = t. + h, j = 1(1)N, 

l J Jo J J *" 1 

and u!?^ 2 is the computed solution at the same point on the mesh 

r - i 2N - - 

| t ± | / where t ± = t^ + h/2, l = 1(1) 2N, 

In the similar way, we can define Zj^ 2 by replacing h by h/2 and N 


by 2N i . e . , 


Z h/2 = m * X I U j /2 ‘ U j /4 


j = 0, . . . , 2N-1 
Now, the computed order of convergence is defined as: 


log Z h - log Z h/2 ] 

Order = 

log 2 

We have taken h = 1/5 and A = 1/16 for finding the computed order 
of convergence and results are shown in Table 3.4. 



Table 3 . 1 


Example 3.1 

h 


-3 

e = 10 

-5 

e = 10 



X p = -1+2E-2 

x = -1+2E-4 

P 


"''max 

1.69395 E -6 

1.69395 E -6 

1/5 

"Smiii 

2.53985 E -12 

2.82659 E -11 


0 

max 

3.37505 E -9 

3.37505 E -9 


0 . 
min 

0.0 

0.0 


*max 

1.04895 E -7 

1.04895 E -7 

1/10 

I . 
mm 

2.46469 E -14 

5.26245 E -13 


0 

max 

3.73001 E -9 

3.73001 E -9 


0 . 
min 

0.0 

o 

• 

o 


"^max 

6.53482 E -9 

6.53482 E -9 

1/20 

*min 

1.54154 E -12 

6.47732 E -11 


0 

max 

3.92125 E -9 

3.92125 E -9 


0 . 
min 

o 

• 

o 

0.0 
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Table 3 . 2 




Example 3 . 2 


h 


-3 

e - 10 

x = 0.4E-2 

P 

-5 

e = 10 

x = 0.4E-4 

P 


*max 

6.79166 E -4 

5.18669 E -5 


^in 

1.55150 E -5 

4.50113 E -7 

1/5 


0 

max 

3.33692 E -4 

4.53325 E -6 


0 . 
mm 

o 

• 

o 

0.0 


"Snax 

3.89329 E -5 

3.17298 E -6 


*min 

3.14849 E -6 

5.34096 E -8 

1/10 


0 

max 

3.32549 E -4 

3.33874 E -6 


0 . 
mm 

o 

• 

o 

0.0 


I max 

3.86764 E -5 

1.36592 E -6 

1/20 

^in 

2.31183 E -6 

2.15286 E -9 

0 

max 

3.32249 E -4 

3.33199 E -5 


0 . 
mm 

0.0 

0.0 
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Table 3.3 




Example 3 . 3 



h 


c = 10“ 3 

x = 0.5E-2 

P 

-5 

e = 10 

x = 0.5E 
P 

-4 


^max 

5.65022 E -4 

5.70341 E 

-5 

1/5 

*min 

1.12327 E -4 

1.32947 E 

-7 


0 

max 

0 . 
mm 

4.97643 E -4 

0.0 

4.90557 E 

0.0 

-6 


"^max 

5.61162 E -4 

5.95588 E 

-6 

1/10 

I . 
mm 

4.97209 E -5 

1.29740 E 

-7 

0 

max 

0 . 
mm 

4.97732 E -4 

0.0 

4.99458 E 

0.0 

-6 


■^max 

5.59555 E -4 

5.93288 E 

-6 

1/20 

■^min 

2.45255 E -5 

4.01312 E 

-7 


0 

max 

0 . 
mm 

4.97737 E -4 

0.0 

4.99925 E 

0.0 

-6 



Table 3.4 


Example 3.1 


Example 3 . 2 


Example 3 . 2 


h 

h/2 

z h 


Order 

1/5 

1/10 

1.49 

E -14 


1/10 

1/20 

1.20 

E -15 

3.63 

1/5 

1/10 

4.7050 

E -4 


1/10 

1/20 

2.5100 

E -5 

4.22 


1/5 1/10 7.0721 E -5 

1/10 1/20 3.6132 E -6 


4.29 



CHAPTER IV 


Variable Mesh Difference Scheme For Linear 
Singularly Perturbed Boundary Value Problems 

4.1 INTRODUCTION : 

Singularly perturbed, boundary value problems occur m many 
areas of engineering and applied mathematics. In many practical 
problems the coefficient of the second derivative is small as 
compared to the coefficient of the first derivative. Examples of 
these are heat transport problems with large Peclet numbers, 
Navier stokes flows with large Reynold numbers etcetra. Because of 
the presence of 'Boundary layers', difficulties are experienced m 
solving problems of above type using numerical methods with 
uniform mesh. In order to get good approximation a fine mesh is 
required in the boundary layer region. 

In this chapter, we first consider the following linear two 
point boundary value problems without first derivative term 

c y" = q(x) y + r(x) (4.1) 

y(a) = a Q , y (b) = (4.2) 

where q(x) > 0 and r(x) are sufficiently smooth functions, ct Q and 
are constants, 0 < e « 1. 

Under these assumptions, problem (4.1) -(4.2) has a unique solution 
and exhibits boundary layers at both ends of the interval (Ascher 
et. al 1988) . 
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(i) S (x) coincides with a polynomial of degree three on each 

[ x i ,x i+l ] ' i= °( 1 ) N - 1 

(ii) S(x) s C 2 [a,b] 

(iii) S(x i ) =y(x i ), i=0(l)N 


The cubic spline can be given in the form (Ahlberg et. al.1967). 

3 .. . 3 


(x.^,- x) 


h?M 


r 11. n n f X. , „ “X \ 

[r<v --P](^ET-) 


S(x) = — 

6h7 M i 

i ' 6h . 

l 

- M. L „ + 
1+1 

+ 

[ y( x i+l> 

- t i+1 ] 

(■?) 



X. 3 
1 

x * x i+i 

where 

- S" ( Xi ) , 

i=0 (1)N 



(4.5) 


It's first derivative is given by 
S' (x) = - M. 2h 


( x i + i- x ) 2 , „ < x - x i> 2 , y( x i + i> - y< x i> w M x „ 

- + M . . „ + r “ n. 


i+l 2h^ 


h i 


x i+l' 1 = OfiJN" 1 


(4.6) 


and the second derivative is given by 


X i+1 X 


S"( X) = + M i+1 -^i 


(4.7) 


For the one sided limit of the first derivative, from (4.6) we 
have 


h i-i h i-i y(x i )-y( x 1 . 1 ) 

S' (X.-) = -±-^= M. . + — M. + 

1 6 l — l o 1 


and 


S'( Xi +) = 


h i h i 

r 1 M i “ -r M i + i + 


i-i 

y(x i+1 ) - y(x.) 


h i 


(4.8) 


(4.9) 
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We have derived a second order convergent variable difference 
scheme using cubic spline for the solution of (4.1) — (4.2) • Then by 
taking second order approximations of first derivative terms, the 
method is generalised for the general linear singularly perturbed 
two point boundary value problem of the form: 

c y"- p(x)y' + q(x) y + r(x) (4.3) 

y (a) = a 0 , y (b) = (4.4) 

where p(x), q(x) and r(x) are sufficiently smooth functions. Under 
these assumptions, problem (4. 3) -(4.4) has a unique solution and 
exhibits boundary layers at one or both ends of the interval 
depending upon the properties of p(x),(Ascher et. al 1988). 

The main idea is to use the 'condition of continuity' of 
cubic spline as a discretization for (4.1) and (4.3). The use of 
cubic spline for the solution of regular linear two point boundary 
value problems without y' term was proposed by (Bickely 1968) . Our 
scheme reduces to that of Bickley's for the corresponding regular 
problem on uniform mesh. 

Also, piecewise C [a,b] approximate solution for off -nodal points 
can be obtained without any further computation. 

4.2 DIFFERENCE SCHEME FOR PROBLEMS WITHOUT Y' TERM : 


Let x« = a 


x . 

l 


= a + 


i-1 

E 1 

k=0 


V 


i=l(l)N, h, = 


x k+l x k' X N 


= b 


Given the values y (x Q ) ,y (x^) , . . . ,y (x N ) of a function y(x) at the 
nodal points x Q , x^, . . . , x N , there exists an interpolating cubic 
spline S(x) with the following properties : 
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(i) S (x) coincides with a polynomial of degree three on each 
[x i' x i+l ] ' ia= 0(l)N-l 

(ii) S (x) e C 2 [a / b] 

(iii) s ( x i) =Y( x i )/ i=0 (1)N 


The cubic spline can be given in the form (Ahlberg et. al.1967). 


_ , v l*rJl 

s < x > “ — ShT 


(x i+l" x)3 (x " x i> 3 

" - M i + ~ 6h. M i +1 + 


h?M. 




[ 


h?M. 


X i * X * X i+1' 1 = 0(1)1X-1 
where M. = S" (x^) , i=0(l)N 

It's first derivative is given by 


+ I *< x i+l> “ 


(4.5) 


(x i+i’ x) , „ ( x ~ x i> . y( x i+i) - M i+r M i 


S/ < X > = " M i -2*- + M i 


i+1 2h i 


h i 


h 


X i " X " x i+i' i = 0 ( 1 ) N “ 1 (4.6) 


and the second derivative is given by 


x i+l" x 


S" (x) = K ± 


+ M. 


x- x i 


i+l h^ 


(4.7) 


For the one sided limit of the first derivative, from (4.6) we 
have 


^i-l ^i-i 

S' (X.-) = -p M.^ + -3— M. + 


yfx.J-ytx.^) (4 . 8) 


1-1 


and 


S'(X.+) = 


h i h i 

- 3 1 M i - ~r M i + 1 


y(x i+1 ) - y(x i ) 


h i 


(4.9) 
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By (4.5) and (4.7), the functions S(x) and S" (x) are continuous 
on [a, b] and for s (x) to be continuous at the interior node x., 

i' 

we have from (4.8)-(4.9), the following condition known as the 
'condition of continuity' 


h i-i h i + h i-i h i 

M. . + — M. + ~ 

6 l-l 3 16 


_ y( x i+ 1 )-y(Xi) yfv-ye’W 

' i+1 * 


h i-l 


i=l ( 1) N-l 


(4.10) 


The 'condition of continuity' ensures the continuity of the first 
derivatives of the spline S(x) at the interior nodes. 

Multiplying and dividing the terms containing h._ 1 by h. and 
h i 

putting a ^ ^ we get 


— jj. 

60^ l-l 


h i 1 h i 

+ "f* 1 + £.> M i + er M i + 1 


y( x i+i) - y( x i) 


°i [ ( y ( x i>- Y(x i-l )] 


(4.11) 


Substituting the following in Eg. (4. 11) 


e M . = q(Xj) y (Xj) + r(x ), j = i,i±l 


we have 


h i \ 1 

65. ^i-iy( x i-i) + r i-i>3 + — (1 + 5 :.)[(q 1 y(x 1 ) + *i)] 


+ T c<J i +1 y< x i + i) + r i +1 5 = 


e[ y (x i+ l)- y(Xi)] eo- i [y(x i )-y(x l _ 1 ) ] 


(4.12) 
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Neglecting the truncation error, we get the following tridiagonal 
system for approximating y lf y 2 / * • ‘^N-l noc * a l points 

x l' x 2' * ’ * ,X N-1 


h i 

e cr . + -= — q. 
l 6<x. ^l-l 


+ 


h‘(l+<r.) 

c( 1 + <7 1 ) + -y q 


h2 i 1 

*i+i = ‘ h i 

[ r i-i + 

(Iter.) r i+1 

- 

6 q i+l 

6<r i 

3 a i 6 



i=l(l)N-l 


(4.13) 


with y 0 - 0 £ 0 , y N = « x 

where = qfxp , r.^ = r(x i ) , i = 0(1)N, 0 ^= h i /h i _ 1 , i=l(l)N-l 


Remark 4.1 

For e = 1 (regular problem) and cr^s 1 (uniform mesh), we have 


r h 2 


r „ , 4h 2 i 

r , ± h 2 1 

[-1 + s 9i_i 

yi-i + 

2 + 6 ■ «i >i + 

r l+ 6 «i + ij 


2 

= - 4 - [ r i-l + 4 r i + r 1 + 1 J. i=KDN-l (4.14) 

y 0 = cc 0 and y x = « x 


is Bickley's scheme for the regular problem 

Y" = q(x) y + r (x) 
y(a) = a 0 , y(b) = 


which 


(4.15) 

(4.16) 
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4.3 DIFFERENCE SCHEME FOR PROBLEMS WITH Y # TERM : 

In this section, we have generalized the scheme derived m 
previous section by taking some second order approximation of the 
first derivative for the problem (4.3)-(4.4). 

Taking the usual Taylor series expansion for y around x^ and 
neglecting the terms containing third and higher order terms, we 


get the following approximations for Y^ +1 and y ^ : 


h? 


i+1 


“ Yi + + -f y'i 


(4.17) 


h? , 

yi-i “ *i - h i-in + y r 


(4.18) 


2 

Multiplying (4.18) by and subtracting it from (4.17), we get 
the following approximation for y^. 


yi “ h i( 1 4^7) [ y i + i + ~ 1)y i - <r i y i-i ] 


(4.19) 


Multiplying (4.18) by and adding it to (4.17), 

following approximation for y" 


y'. ' - 
J l 


2<r i 


h i (1+<r i) 

Also we have 


[ y i+ i - + + Vi-i . 


we get the 


(4.20) 


y 7 . 7 
■1 y l 


yi+i “ y i + h i y i' 

y i-i “ y i - h i-: 

Using the expressions for y^ and 
respectively and putting them in 
approximation for y^ +1 , 


(4.21) 

(4.22) 

y^ 7 from (4.19) and (4.2 0) 
(4.21) we get the following 
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y i+i “ HrfelT [ (2<vi)y i+1 - (VDS + *Iy l _ x ] 

i=l(l)N“l (4.23) 

Similarly, using the expressions for and yV from (4.19) and 

(4.20) respectively and putting them in (4.22) we get the 
following approximation for , 

*t-i “ h^ti+^T [ " y i+i + (CT i + 1)2 + 2 ) yi-J 

i=l(l)N-l (4.24) 

Substituting the following in Eq.(4.10) 

S «j = P(Xj)yj+ qCx.-Jy.. + r(x..), j = i,i±i 

and using Eqs. (4.19), (4.23) and (4.24) for the first order 

derivatives, we get the following system which gives the 
approximations for y , y 2 , ..., y N-1 . 


h? 

e <r + 
l 6cr 


h.(ai+2) 


Vi 


. 2 
h <r 
l l 


. q i-l 6(1+0-.) p i-l “ 3 P 1 + 6 ( l+o- 1 ) p i+l 


Vi 


h?(l+cr.) h. (1+cr.) h. (<r?-l) h.(<r.+l) 

e(l+«r i )+ < Ji + g FT - p i-l 3a\ p i 6 p i+l 

111 


h? 


h. 


h. 


- c + if- q 


h. (2<r.+l) 

i l _ 1 ^ l v i ^ 

“ q i+l “ 6 (cr^+1) cr^ p i-l + 3 o* i p i 6 (o^+l) p i+l 


i+1 


*1 


r i-i (1+cr i ) x r i+i 

+ ± r . + 

6c r. 3cr. l 6 


i=l (1) N-l 


(4.25) 


with y 0 = « 0 f y N = «! 
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where p i = p(x i ), f i - f(x.), i = 0(1)N, cr.= hj/h^, i=l(l)N 

For selecting the mesh points, we follow the procedure given by 
Jain et. al(1984). 

4.4 MESH SELECTION PROCEDURE : 


Let N be the number of mesh points in the interval [a,b] and 
(=h^/h^_ 1 ) be the mesh ratio (as defined earlier). For simplicity 
we take cr = cr ( a constant) V i. 


We have 


b - a = x N - x Q 


= ( X N _X N- 1 ) + < X N-l' X N-2 ) + "- + + ( x r x 0> 

= cr N ^h.. + cr N 2 h Q + . . . + crh^ + h Q 


,N-1 N-2 i \ u 

= (cr + cr + . . . + cr + l)hg 


which implies 


h o = 


(b-a) (o-l) 
(cr N -l) 


{ (b-a) (l-<r) 

(l-o* ) 


if cr > 1 


if cr < 1 


(4.26) 


Therefore, given the values of N and cr, we can choose h^ from 
(4.26) and subsequent h^'s can be obtained as = crh^^, i=l(l)N-l 

CASE I 

If the boundary layer occurs at the left end then we choose 
cr > l. This ensures more mesh points near the left end of the 
interval . 
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CASE II 

If the boundary layer occurs at the right end then we 
choose cr < 1. This ensures more mesh points near the right end of 
the interval. 

CASE III 

If there are boundary layers at both ends then we consider 
first half of the interval [a # (a+b)/2] with o > 1, then we take its 
mirror image in other half interval [(a+b)/2,b]. This ensures more 
mesh points at both ends. 

CASE IV 

If there is boundary layer at the center then we consider 
first half of the interval [a,(a+b)/2] with cr < 1 , then we take 
its mirror image in other half interval [(a+b)/2,b]. This ensures 
more mesh points at the center. 

4.5 ERROR ANALYSIS : 

Taking the usual Taylor series expansion for y around we 

get the following expressions for y^ +1 and y 
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Multiplying (4.28) by cn, and subtracting it from ( 4 . 27 ), we get the 
following expression for y^ 


h? 


y'i = y,(x i> + 6 ti%) [ + £ y'"(^ u ) ' 


(4.29) 


Multiplying (4.28) by and subtracting it from (4.27), we get the 
following expression for y^' 

y v =y"(x i ) + -gi [y'"(e^ u ) - -J- y"'(Z l 2 l) ) ] (4.30) 


Also, we have 


y 7 * . i - y- + h.yV 
1 l+l J i i J i 


h? 


(i) 


y'" ($;■ ) 


(4.31) 


? 3 l> “ X i+1 


Putting the values of y^ and y^ ' from (4.29) and (4.3 0) in 
(4.31), we have 

9 9 

2 + cr. . 


n? 


n? 


y^i = y'( x ^i) + 


ii . . ii • r r • w * n 

_i y"'(£ (1 >) + -L- i 

2 j * »3 7 6 \_( 1 + cr^ j 


(-? — - )y'" (Cj 11 ) ] 

Vf(cr.+lW J 


+ 


i( cr i + l) 


y'"(€i u ) 


(4.32) 


Similarly, we have 

y l-i = n - h i-i y i 


" + y"' (€ 4 °) 


(4.33) 


( i) 


X i-1< ^4" < X i+! 
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Putting the values of y^ and yV from (4.29) and (4.30) in 


(4.33) , we have 


h ? 


h? 


y i-i = y ' (x i-i> + ~2 + 4-[(-57W i )) y '"< >) 


2a*. +1 


( 45U . tj l n _ 

-2 y #// (C2 i} ) 1 

0*7 (0* i + 1 ) * 1 J 


(4.34) 


Expressions (4 . 29) , (4 . 32) and (4.34) show that y:, y^ +1 and Y'^ ± 
are second order approximations to y' (x ) ,y' (x. +1 ) and y' (x^_ 1 ) 
respectively. 

Putting the tridiagonal system (4.25) in the matrix form, we 
have 


MY = R (4.35) 

where M = (m. ), 1 as i,j < N-l is a tridiagonal matrix with, 

-*- / J 

m^ i+1 = coefficient of Y^ +1 m (4.25) i=l(l)N-2 

m. . = coefficient of y. m (4.25) i=l(l)N-l 

1,1 i 

m . = coefficient of y. n in (4.25) i=2(l)N-2 

and R =(r^), 1 ^ i < N-l is a column vector with 


h 2 


2(l+cr i ) 


n i T 1 z '- LTU i ' 1 

r i - “ T- r i-l + trT r i + r i+l J’ 


and Y = ( y 1( y 2 , . . . ) 


tr 


i=l (1) N-l 


Also we have 

MY a - T(h) = R (4.36) 

where Y^= (y (x^) ,y (X£) , • • • tY ( x fj-l) ) ^ denotes the actual solution 
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and T (h) = [ T^) ,T(h 2 ) , . . . ^(h^)! 


tr 


E0 ’i 4 (iv) 

where T^) = ($.) 


x i-l< «i< x i 


l+l 


is the truncation error 

From (4.35) and (4.36), we have 

M(Y-Y ) = T(h) 

Thus, the error equation is 

ME = T(h) 


where E = Y - Y a 

A 


(4.37) 


(4.38) 


(4.39) 


Also, it can be seen that for h^ satisfying following conditions 
for different cases, the tridiagonal matrix M is irreducibly 
diagonally dominant. 


Case I 


p(x) > 0, q(x) < 0 


For m. ._ < 0, we must have 

1,1 l 


A., + h7A 0 - h A, - h.A, + h.A,. < 0 

1 12 13 14 ID 


(4.40) 


where 


A 1 = e (Ti 


A 2 6<r i q i-l 


A 3 6(l+<r i ) P i-l 


cr . +2 
l 


rr 



A 5 6(l+<r i ) P i+l 
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Since A x > 0, it is enough to have 

h i A 2 “ h i A 3 ” h i A 4 + h i A 5 < 0 (4.41) 

Since A 2 < 0, this implies 

-i-( A 3 + A 4 - A 5 ) < h. (4.42) 

Let H x = -jj- ( a 3 + A 4 - A 5 ) 

Th6n H 1 = (iX - ^ - q — 1 [ (J x +2) Pi-l + - ^Pl + l] < 4 ‘ 43 > 

Therefore, for m. . - < 0 , we must have 

i , i — i 

H 1 < h. (4.44) 



(4.46) 
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Since B 2 < 


Let H 2 = 


Then H 2 = 


Therefore, 


For 


where 


Since C 1 > 


Since C 2 < 


0, this implies 

h. < - 4- f 

1 B , l 


B 3 - B 4 - B 5 


(4.47) 


B , ( 


B 2 ( B 3 ” B 4 B 5 


) 


[Pi-1 + 2 «V1) p ± - <T iPi+1 ] 


(4.48) 


for m. . > 0, we must have 


h i < H 2 


(4.49) 


m . - < 0, we must have 

JL g JL i JL 


- C 1 + h i C 2 * h i C 3 + h i C 4 + h i C 5 < 0 


C ! - C 


= 


2 6 ^ 1+1 


(4.50) 


C 3 6(l+or i )cr i P i-1 


C . = 


4 3 cr *1 


C 5 = 6(l+<r i ) P i+l 

0, it is enough to have 

h?C„ - h.C 0 + h.C. + h.C K < 0 
12 i3 14 15 

0, this implies 

\ ( C 3 ' C 4 - C 5 > < h i 


2<r i +1 


(4.51) 


(4.52) 



Let H 3 = g- ( C 3 - C 4 - c 5 ) 


Then H 3 = — 


i (i'^ i ')q i+1 [ P i-l + 2{1+a i ] Pi + 2 < 1+ff i> P i+ i] ( 4 - 53 ) 


Therefore, for m i < 0 , we must have 


H 3 < h i 


(4.54) 


Combining the conditions (4.44) , (4.49) and (4.54), we have 

max(H 1 ,H 3 ) < < H 2 (4.55) 

Similarly, for other cases, we have the following conditions 
Case II p(x) > 0, q(x) > 0 

H 1 < h x < H 2 (4.56) 

Where H„ = - ~ — \ p. , + (cr.-l) p. - 2ar p. ,1 


Case III 


where 


P(X) 

> 0, q (x) 

H 1 < 

h i < H 2 

1 


*1- 

-! L Pi -! + 

cr. 

1 

r °V 2 


«i+i 

L cn+l p i- 

POO 

< 0, q (x) 

,k 2 ) < 

h i < H 3 

•H 

b 

VO 

1 

r "i! 2 p 

*i-l 

L<r i ((r i + 1 ) p 

1 

J 1 


(4.57) 


2 2<r^+l -i 

~ q i-l [o' i (o' i +l) P i-l ~ ?i cr i +1 Pl+1 J 


" 3£r i T °i _1 °i +1 cr i +1 I 

H 3 ~ (1+cr.^q.^ 3<r^ p i " 6 p x+l + 6 a p i-lJ 


and 
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Case IV 


p(x) < 0, q(x) > 0 
H x < h 1 < min (H 2 , H 3 ) 


Where 

H i 

= 4l 

cr.p. . , + 
1*1+1 

(1 * 

a i )p i- Pi-1 


H 2 

<T . 

1 

r °" i+2 p 


cr? 

2<r i p i" IhI.1 
1 


H 

1 

•H 

tJ 1 

L ^i+i p 

i-l + 


and 


Case V 


where 


Case VI 


and 


i r i , 2£r i +1 l 

[ i+l L c V a 'i +1) Pi_1 " Pi ~ °‘i +1 P i +1 J 


H. = 


p(x) s 0, q(x) > 0 
h^ < min (H.^ H 2 ) 


H i = 


o' i / 6 e 


and 


6 e 


«i-i 


H„ = 


*i-l 


p(x) < 0, q(x) s 0 

< min (H 1# H 2 , H 3 ) 


ecr 


where H„ = 


1 <r ( 1+cr^) P i+1 

6cr^c (o^+l) 


HL = 


2 2(o‘ i +l)p.+cr i (2cr i +l)p i+1 -p i „ 1 


HL = 


6cr^c 


3 ^iPi+i" 2 (cr i" 1)P i" P i-l 
Case VII p(x) > 0, q(x) 3 0 


(4.58) 


(4.59) 


(4.60) 


h^ < min (H^, H 2 , H 3 ) 


(4.61) 
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Case IV 


p(x) < 0, q(x) > 0 
H 1 < h i < min (H 2 , H 3 ) 


Where 

H i 

= 4L 

cr . p . , „ + 
1*1+1 

(1 - (T.jp.- p._ 1 


H 2 

cr . 
i 

r p 

J_ O/r v\ - 


«i-i 

L P 

l-l l*i l+cr. J 

l 


and 


Case V 


i r i , 2 <r i +1 i 

[ i+1 L cr i< cr i +1 ) P i _1 " Zp i” O'i+l P i+lJ 


H„ = 


p(X) a 0, q(x) > 0 

h^ < min (H 1# H 2 ) 


where 


Case VI 


where 


H i = 


cr i / 6e 


and 


<2h- 


HL = 


l-l 


p(x) < 0, q(x) s 0 


hu < min (E^, H 2 , H 3 ) 

ccr i 

H 1 = oTTl+oT) P i+1 


= 


6o* i e(cr i +l) 


2 2(cr i+ l)p. + cr i (2(r i +l)p i+1 -p i . 1 


and 


Case VII 


= 


6cr^e 


3 °‘i p i+l “ 2 (<r i _ 1 )p i" p i-l 

p(x) > 0, q(x) s o 


6c 


*i-l 


(4.58) 


(4.59) 


(4.60) 


h^ < min (H 1# H 2 , H 3 ) 


(4.61) 



Case IV 


p(x) < 0, q(x) > 0 


and 


H x < h^ < min (H 2 , H 3 ) 


Where H- = 


= 


2 q 


= 


|.[ *1*1+1 + <* - Pi-J 

r 0‘ i +2 ■ a 2 , 

r^l ft+i p i-i + 2<r iPi- T+s.Pi+1 J 


1 1 2a i + l 1 

[ i+l L cr i( tr i +1 > Pi_1 ~ P i” o-i+l p i+lJ 


Case V 


p(x) S 0, q(x) > 0 


< min (H 1# H 2 ) 


where 


H i - 


cr. / 6e 


and 




= 


1-1 


6e 




l-l 


Case VI 


where 


and 


p(x) < 0, q(x) a o 


h^ < min (H 1 , H 2 , H 3 ) 


eo i 

H. = - P. 


1 cr i (l+cr i ) *i+l 


H„ = 


6cr^c (cr^+1) 


2 2( £ r.+l) P .+<T.(2a i +l)p. +1 -p i _ 1 


H. = 


6<r. e 


3 <r i p i+r 2 (<r i _ 1)p i _p i-i 


Case VII p (x) > 0, q(x) = 0 


(4.58) 


(4.59) 


(4.60) 


h^ < min (H 1 , H 2 , H 3 ) 


(4.61) 
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N-l 

e. = £ »T^T(h) 
1=1 


j=l (1) N-l (4.67) 


and therefore 


e . ^ — 

J 1 i-2 


Ch 


h l i°i 

o o 


j=l (1) N-l (4.68) 


where C is constant independent of h = max h. 

l=sisN-l 1 

Therefore, 


E I = 0 (h*) 


h = max h. 
l<i<N-l 1 


(4.69) 


4.6 NUMERICAL EXAMPLES : 


Example 4 . 1 (Doolan et. al. 1980) 

e Y" = y + cos 2 ir x + 2n 2 cos 2 ttx 

y(°) = y(i) = o 


The exact solution is given by 


, v = exp(-(l+x)//g) + exp(-x//c) 
Y ' ' 1 + exp (-l/v'e) 


2 

COS 7TX 


Since p(x) =0 and q(x)=l > 0, the boundary layer exists at both 
ends, so we choose mesh according to case III of the section 4.4. 


Example 4.2 (Jain et. al. 1984) 
ey" = y' 

y (0) = 1, y(l) = 0 


The exact solution is given by 


y(x) 


1-exp (-(1-x) /e) 
l-exp(-l/e) 
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Since p(x) ■ 1 > 0, the boundary layer exists near the right end, 
so we choose mesh according to case II of section 4.4. 

Example 4.3 (Reinhardt 1980) 

c y" = - y' + l + 2x 

y(0)=0 , y(l)-l 

The exact solution is given by 

y(x) = x(x+l-2e) + (2e-l) l~exp(~x/e) 

' 1-exp (-1/e) 

Since p(x) a -1 < 0, the boundary layer exists near left end, so we 
choose mesh according to case I of section 4.4. 

4.7 DISCUSSION : 

The second order variable mesh difference scheme, obtained 
from the 'condition of continuity' of the cubic spline, can 
be applied to any linear singularly perturbed boundary value 
problem, provided that there is some prior information about the 
location of boundary layer, so that mesh can be chosen 
accordingly. 

We have implemented our method on three examples. Maximum 
error at the nodal points, max | y(x^) - y^ | , are tabulated in 

Tables 4. 1-4. 3 for different values of the parameters e, N and 

cr(for simplicity we have taken cr^= cr V i ) . 

It is observed that errors are less for variable mesh (cr * 1) 

compared to that of uniform mesh (or = 1) , which shows the need of 

the variable mesh for the problem with boundary layer region where 
the solution changes rapidly. The uniform mesh gives good results 
when mesh is taken very fine through out the interval, which is 
quite demanding on the computer time. Also solution at any 
off-nodal point can be obtained from the cubic spline expression. 
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Table 4.1 


Max. absolute error (Example 4.1) 







a 






N 

1.00 


1.05 


1.15 




60 

1.87947 

E +1 

8.29368 

E -2 

3.81489 

E -3 

e = 

10 5 






. 




120 

9.54587 

E -2 

3.63303 

E -3 

6.63985 

E -4 


-8 

100 

2.67670 

E -1 

2.56806 

E -1 

9.78686 

E -3 

e = 

10 

200 

2.66837 

E -1 

5.97463 

E -2 

6.43276 

E -4 


„ -10 

150 

2.67942 

E -1 

2.60797 

E -1 

2.02680 

E -3 

e = 

10 

300 

2.67924 

E -1 

4.69272 

E -2 

6.41468 

E -4 


„ -12 

200 

2.67949 

E -1 

2.61329 

E -1 

9.44750 

E -4 

e = 

10 

400 

2.67948 

E -1 

3.66378 

E -2 

6.40522 

E -4 
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Table 4.2 


Max. absolute error (Example 4.2) 







0* 





N 

1.00 



0.95 


0.85 


~5 

60 

6.97322 



2.37996 


5.14372 

E “2 

e = 10 

120 

1.91001 



8.08277 

E -1 

3.20425 

E -3 

e = 10' 8 

120 

1.73610 

E 

+3 

1.99962 


4.14263 

E -3 


240 

4.34024 

E 

+2 

1.98215 


3.19924 

E -3 

e = IQ' 10 

250 

4.00000 

E 

+4 

1.96678 


3.20501 

E -3 


500 

1.00000 

E 

+4 

3.38173 

E -4 

3.20504 

E -3 

. “12 

300 

2.77777 

E 

+6 

1.96644 


3.22245 

E “3 

O 

II 

H 

O 

600 

6.94444 

E 

+5 

7.33741 

E -5 

3.22245 

E -3 




400 


3.12550 E +6 
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e 4.3 

) 

<r 

1.05 1.15 

2.35713 6.09777 E -3 

1.48059 2.36813 E -3 

2.01874 2.83734 E -2 

1.97420 E +2 2.36227 E -3 

1.98307 6.39184 E -2 

1.96816 2.36184 E -3 

1.97073 4.00547 E -3 


1.96800 


2.36315 E -3 



CHAPTER V 


Third Order Variable Mesh Methods For Semi-Linear 
Singularly Perturbed Boundary Value Problems 

5.1 INTRODUCTION : 

Consider the following class of semi-linear singularly 
perturbed boundary value problems 

ey" = f(x,y) a<x<b, 0 < e « 1 (5.1) 

y (a) = a, y (b) = £ (5.2) 

In many fields of application, including geophysical fluid 
dynamics (Carrier , 1970) , these types of problems arise which 
exhibit boundary layer of small thickness at both ends. Difference 
methods with uniform mesh are not suitable to these type of 
problems as a fine mesh is required in boundary layer region and 
comparatively much coarser mesh elsewhere. 

It is known that, under suitable assumptions, one expects the 
solution y(x) of (5.1) -(5. 2) to behave qualitatively as the 
solution of the following problem (Howes, 1980) . 

ey" = p(x)y + q(x) p(x) * p > 0 (5.3) 

y (a) = a , y(b) = £ (5.4) 

where p(x) and q(x) are sufficiently smooth functions. 
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In this chapter, a two parameter family of methods based on 
cubic spline approximation, which gives third order approximation 
to the solution of (5.3) — (5.4), is presented. The methods are then 
used for solving (5.1) -(5.2) by quasi-linearization technique. For 
approximating solution at the nodal points, we have given a family 
of third order variable mesh methods. In addition to the parameter 
o*^, one more parameter X is introduced which is exploited to yield 
higher order approximations of the solution at off-nodal points. 
These third order variable mesh methods reduce to a one parameter 
family of fourth order methods when the mesh ratio is taken to be 
equal to unity i.e for the uniform mesh. In addition, when the 
parameter X is equal to unity then we have well known Numerov's 
method for the corresponding regular problem. Also, the restriction 
on mesh size for convergence can be relaxed by suitably choosing 
the value of the parameter X. 

Earlier, Jain et al. (1984) have given a one parameter family 
of third order methods for the solution of y" = f(x,y,y') where f 
contains the small parameter c implicitly. These methods reduces to 
one particular method namely Numerov's method for uniform mesh. 
Their methods become particular cases of our methods when the 
parameter X = 1 for the problem of the form (5.3)-(5.4). 

In order to get family of higher order methods instead of one 
particular second order method, as derived in the previous chapter, 
we proceed in the following way. 
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5.2 DERIVATION OF THIRD ORDER METHODS : 


Let x Q = a, x 4 = a h k = x k+l " x k' X N = b ‘ 

Consider the Ahlberg's cubic spline, considered in the previous 
chapter . 


<x 1+1 - X) 

S(X) gjr 


M i + 


(x - x.) 


h?M. 


- r 1A i rA i 1 ( x i+i” x 1 

6h— M i + 1 + L Y(X i> - — J l fa— ] 


h 2 M. 


r , , “i"i+i i r x_x i i 

+ [ y( X i+l) " — 6 J 1 — h— J 


x i s x * x i+l' 1 


= 0 (l)N-l 


(5.5) 


where = S" (x^) , i=0(l)N 


Now, consider the following difference scheme for computing the 
approximate solution ^ 1 *^ 2 ' * * * '^N-l the nodal points 

X 1 ,X 2 ' * * * /X N-1* 

- c * 01 * 1-1 - c a n y i " e y i+i 

h l[ B oin~x + B n y i' + B 2i y i;x] - 0 (5 - 6) 

i = 1(1) N-l 


+ e 


with 

y 0 = y(a) = a , y N = y(b) = |3 


(5.7) 


I < x s 1 ' 


where 
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The coefficient a's and B's are determined by using Taylor's series 
expansion and equating the coefficients of h^y (r) ( X;L ) , r = 0(1)4 to 


zero. 


Thus, we have the following set of five equations: 


“ a 0i + a li - 1 = 0 


(5.8) 


-1 = 0 


(5.9) 


ZJ a 0i - I + < B 0i + B li + B 2i > = 0 


(5.10) 


6o .3 a 0i " 6 " cr B 0i + * B 2i ~ 0 


(5.11) 


24<r? 301 
1 


1 ^ 

24 + ^2 B 0i + 2 B 2i 


(5.12) 


where cr . = h . /h . . 

l l l-l 


From (5 . 8) - (5. 12) , we have 


a oi = 


(5.13) 


a ii = 1 + cr. 


(5.14) 


^(cr.- 1) (1-2X) + 1 


12 A cr i 


(5.15) 


(erf + 1) ( 2A-1) + 2A(3X-l)cr i ( cr.j+1) + (2A-1) 


2 2 
12 a erf 


(5.16) 
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and 


B 


<t ± + (2A-1) (<r i ~ 1) 


2i 


12 


> 2 rr 2 

X <T i 


(5.17) 


For c yV and e y^+ x# w ® take the following approximations 

e Y-' = P i Y i + q ± (5.18) 

and e y^ x = P i±x S(x. ±x ) + q. ±x (5.19) 

where p. ±A = p(x. ±x ), p. = p(x.) , q i±x = q(x. ±x ) ,q x = q(x.) , 
and S(x... ) are obtained by putting x = x +Ah and x = x -Ah. - 

l-A 11 1 1 _ 1 

in S(x) given by 


(Xi - x) 3 _ (X - X .) 3 _ 

“ M. + 

1 - 6h i i+1 


s < x > sir M i + 



x i+r x 1 

L y i 6 J 

h. 1 

L 1 J 


h?M. 


r “i"i+i i ( x “ x i i 

L y i+i - — J 1 — h— J 


x i ^ x ^ x i+1# i = 0 ( 1) N-l 


(5.20) 


where 


e Mj = p(x.)y . + q(Xj) 


3 = i,i+l (5.21) 


Putting the values of c y^' and e along with values 
of a' s and B's, we get the following tridiagonal system which gives 
the approximate solution y 1 ,y 2 * * * * ,y N-l at the nodal points 

x l' x 2 ' • * * x n-1 
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[- E<r i + h i B 0i p i-X ( X + p i-i a l( X ' h i-i>)] Vi-x 
+ [ e(l + <r.) + h\ p. B ti + h*(l-X) (b 0 . p._ x + B 2 . p. +A ) 

+ h i P i( B 0i p i-x a 2< X ' h i-l> + B 2i P i+X a 2 (X ' h l>) ] y i 

+ [- £ + h i B 2i P i+X ( X + P i+1 a i( X ' h i>)]y i+ l 


= - ^ {a 1 (^h._ 1 ) B oi p._ x < J i _ 1 + a l( X,h.) B 2 . p i+x q i+1 

+ a 2 (X,h i _ 1 ) B oi p._ x q. + a 2 (X,h.) B 2i p. +x q. 

+ B Oi q i-A + B 2i q i+l} " h i B li q i 


i = 1(1) N-l (5.22) 

y 0 = y(a) = a , y N = y(b) = p, (5.23) 

where a x (X,h) = h — * / a 2 (X,h) = h ^ - ^ -[(1-A) 2 ~l] . 

Remark. 5 . 1 

For e = 1 (regular problem), <r^ = 1 Vi (uniform mesh) and X=l, 
we have from (5.22) 


(-1 + 12 Pi-i) Vi-1 


2 + 


10h‘ 

12 


+ 


- 1 + 


12 P i+lJ 


i+1 


h* 

12 


q i-l + 10q i + q 




i = 1 ( 1 ) N-l 


(5.24) 


y 0 = y N = P 
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which is 3. well known fourth order Numerov , s method for the regular 
problem 

Y" = P(x)y + q(x) (5.25) 

5.3 SEMI -LINEAR PROBLEMS : 

For semi-linear problems of the form (5.1)-(5.2), we use the 
quasi-linearization technique. Starting from the initial guess 
u Q = u Q (x), we obtain the sequence of linear boundary value 
problems 

eu n+l = f y (x 'V u n+l + f < x 'V ' V x ' u n )u n ( 5 - 26 ) 

u n+l (a) = u n+l (b) = P (5.27) 

2 

It is known that (Doolan et.al. 1980), if is bounded for 

ay 2 

all x e [a,b] and all real y, then there exists a constant C, 
independent of n and e, s.t. 

2 

U n+1 ( x ) “ u(x) " C u n (x) ” u(x) (5.28) 

where u(x) is the solution of (5.1) — (5.2) . 

We note that each linearized eq. (5.26) can be written in the form 
(5.3) by putting 

p(x) = f y (x,u n ) and q(x) = f(x,u n ) - f y (x-V u n (5.29) 

In order to exclude turning points in (5.26)-(5.27), we assume 
that, for each n £ 0, a constant K n exists such that 

f y (x 'V * K n > 0 ' 


x € [a, b] 


(5.30) 
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5.4 MESH SELECTION PROCEDURE : 


For simplicity we take cr^ = cr (V i) . 

Since there are boundary layers at both ends, we first 
consider first half of the interval [a, (a + b) /2 ] with cr > 1, 
Given the values of N ( an even positive number) and cr, we can 
choose h Q from (4.26) of previous chapter as 


= (b-a) (tr-1) 

0 " 2«T N / 2 -l) 


for cr > l 


(5.31) 


and subsequent h^'s can be obtained as h^ = crh.^, i=l,N/2-l 


Then we take its mirror image in other half of the interval 
i.e.,[(a + b)/2, b] . This ensures more mesh points in the boundary 
layer regions at both ends. 


5.5 ERROR ANALYSIS : 

In matrix-vector form, the tridiagonal system (5.22) can be written 
as 

MY = D (5.32) 


where M = (m^), l s i,j s N-l is a tridiagonal matrix with 


m. . , „ = - e 
i,i+l 


+ h i B 2i P i + A ( X + Pi+1 a l (X 'V) 1 = N_2 


m. 


l# i 


= c (1+ cr . ) + 


* 2 i 


p. B 


li 


h 2 (l-X) 


0i 


Pi-; 


+ B 


21 *l+X 


+ h i P l[ B 0l P i-A a 2 (X ' h i-l> + ®2i P i+A a 2 (X ' h i). 


i = 1, N-l 
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= ' e + h i. B 0i Pi-X 0 + Pi-1 Vl)] 


and D — (d^) , 1 - i - N-l is a column vector with 


i = 2, N-2 


d i = ‘ h i[ a l( X ' h i-l> B 0 i Pi-x q i-l + B 2i p i+x q 


1 + 1 


+ a 2 (X,h._ 1 ) B 0 . p._ x q.^ + a 2 (X,h.) B 2i p i+x X q. 
+ B 0i «i-X + B 2i <*i + X + B l‘Ii] 


and Y = . . . -y ^) 1 • 

Also, we have 

MY a - T (h) = D (5.33) 

where Y a = (y (x^) ,y (x 2 ) , . . . ,y (x N-1 ) ) denotes the actual solution 
and 

T(h) - ^T(h^) , . . . ,T(h N _ 1 ) j is the local truncation error with 
T(h^) given by 

TO^) = — £-5 [(5X-3) (o-?-^) + 5X(2X-1) (erf- <r?)lh?y (v) ( Xi ) + 0(h®) 

3 6 Ocr. *- j 

1 

(5.34) 


From (5.32) and (5.33), we have 

M(Y - Y) = T (h) (5.35) 

A 

Thus, the error equation is 


ME = T (h) 


(5.36) 
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where E = Y a - Y. 

From (5.15), it can be seen that B Qi 's are positive for | < A s 1, 

if o*. satisfies 
x l 

o* i (o*i“ 1) (1-2A.) + 1 > 0 (5.37) 


which implies, 


(2A.-1) <r? - (2A.-1) o* i -1 < 0 


(5.38) 


which is equivalent to 



(5.39) 


Therefore, either 



(5.41) 

Since cr i > 0, second inequality in condition (5.40) is always 
satisfied and second inequality in condition (5.41) does not hold 
always, we must have 



(5.42) 
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Similarly , from (5.17) it can be seen that are positive for 

u X 

1 < X * 1, if <j. satisfies 

2 

<r?+ (2X-1) (<T^-1) > o (5.43) 

which is equivalent to 



Therefore, either 


(5.44) 



(5.47) 

Since cr^ > 0, second inequality in condition (5.46) is always 

satisfied and second inequality in condition (5.47) does not hold 
always, we must have 

cr. > (2X-1) (- | + \ 



Combining conditions (5.43) and (5.48), we have 




108 


For the tndigonal matrix M to be irreducibly diagonally 
we do the f ollwing if m.^ 1+2 sre non- positve. 

for m. . - < 0, we must have, 

1 , i“i 

• 6 "i + h i B 0i Pi-X 0 + Pi-1 a i<*. h i_i)] < 0 


Now, if 

h i B 0i Pi-X l> + Pi-1 a i<*< h i-i)] < 0 

then (5.50) is satisfied 

Substituting the value of a 1 (X,h^_ 1 ) we have 


h? 


h i B oi Pi-x l> + Pi-1 4 ir ^ a 2 -d] < 0 

which is equivalent to 


h? 

h i B oi Pi-x x + T~2 B 0i Pi-x Pi-1 < 0 

6ecTi 


which implies 


6 ccrl < h? p^d - X 2 ) 


Similarly, for m. . < 0, we must have, 

X f It jL 

- c + h l B 2i Pi+X ( X + p i+l a l (X ' h i ) ) < ° 

Now, if 

h l B 2i Pi + X ( A + Pi + 1 a i (X ' h i>) " ° 
then condition (5.55) is satisfied 


dominant , 

(5.50) 

(5.51) 

(5.52) 

(5.53) 

(5.54) 

(5.55) 

(5.56) 
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Substituting the value of a 1 (X,h i ) # we have, 

h? 

h i B 2i P i+X C x + Pi+l -T r x < x2 - 1 )] < 0 (5-57) 


which is equivalent to, 


hf 


h i B 2i P i+X X + 6F - B zi P i+X Pi+l < 0 (5.58) 


which implies, 


6c < h i - x2 > 


(5.59) 


Let £ = min p(x) , then condition (5.54) and condition (5.59) is 
a^x^b 


satisfied if 


6e<r? < h 2 P (1-X 2 ) 


(5.60) 


For m. . > 0, we must have 


i,i 


E (l + a x ) + h 2 Pl B u + h 2 (l-X)(B 0i Pi _ x + B 2i p. +A 

+ h i P i[ B 0i P i-X a 2 (X ' h i-l> + B 2i P i+X a 2 (X ' h i)] > 0 

(5.61) 


It is enough to have, 

h i p i B u + h i^- X > ( B 0i p i-x + B 2i p i + xj 

+ h i P i[ B 0i P i-X a 2 (X ' h i-l> + B 2i P i+X a 2 (X ' h i)] 

which is equivalent to, 

Pi Bli + (1 -X) (b 01 Pi _ x + b 2 . P i+X 


> 0 
(5.62) 
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+ Pi[ B Oi P i-X a 2^' h i-l> + B 2i P i+ X a 2< x ' h i)] > 0 


Substituting the value of a 2 (X,h i _ 1 ) and a 2 (X,h.), we have 


(5.63) 


Pi B 1 . + (l-X) ( b 0 . Pi _ x + B 2 . p. +x ) > 

h 2 ^2 
p i B oi p i-x ~tir (i-i)[i-d-x 2 )] + p.B 2i P . +x -|j (1-X)[1-(1-X 2 )] 


(5.64) 


Let T = max p(x) , then it is equivalent to 
a^x^b 

9 ( Bl . + (l-X) (B 0 . + B 2 .) ' > 

2 h? h? 

9 |B oi ~ 2 (l-X) [1-(1-X 2 ) ] + B 2l (1-X)[1-(1-X 2 )] 


6e<j . 
l 


(5.65) 


which is equivalent to 


6eo j£ l B ii + d-^(B oi + B 2 ) 
? 2 X(1-X) (2-X) (b o1 + <r 2 B 2i 


For diagonal dominance, we must have 




(5.66) 


M i,i I ‘ I m i,i 


i,i-l I - I m i,i+l I > ° 


(5.67) 


Substituting these values, we must have 
C(l+ o-.) + h 2 p. B u + h 2 (l-X)(B 0 . Pi _ x + B 2 . p. + J 



Ill 


+ h l Pi[ B Oi Pi-X a 2< X ' h i-l> + B 2i P i+ X ^(X.hJ] 


" E °i + h i B 0i P i-X l> + Pi-1 V X ' h i-l>] 

" G + h i B 2i P i+X ( X + Pi+1 a l ( X ' h i>) > 0 

Which is equivalent to 


(5.68) 


l i ( Pi-X B 0i + Pi B li + P i+ X B 2i ) > 


hf 


6eo^ 


2 X(1 X J ( B 0i p i P i-A +cr i B 2i p i+l p i+A 

(T . ' / 


+ X ( X - X ) < 2 " X > ( B 0 i PiPi-X + °L B 2 i P i+ lPi + x) ] < 5 - 69 > 


Therefore, it is enough to have 


? ( 


B 0i + B li + B 2i j 


h 2 

i -2 2 

> — ~ ? X(l-X ) 
2ea1 


B 0i + 


4 B 2i ) 

(5.70) 


l. e. , 


2CO-? P f 


B 0i + B li + B 2i 


.)> 


f P 2 X(1-X 2 )(b 0 . + (T 2 


* B 2i 


(5.71) 


which implies 
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2 £<?. 


i £ ( 


B 0i + B li + B 2 


i) 


> h? 


Mi-x 2 )(b 01+ cl B 2i ) 


(5.72) 


Combining conditions (5.60) , (5.66) , and (5.72), we have the 
following condition on h^ for M to be irreducilbly diagonal 
dominant and hence monotone, 


where 


H 1 < hi c H 2 H 3 


(5.73) 


H 1 °i -v 2 


£( 1-0 


ecr i~ 


H 2 = 


P 2 M1-X)(b 0 .+ <r B B 2 ^ 


and 


H 3 = min 


6 ( B i + d-*> < B oi + B 2 .)) a( Bli + B oi+ b 2 i j) 


(2-X) 


(1+X) 


from eq. (5.36) 


E = M” 1 T (h) 


(5.74) 


N—1 f 

Let s . = Y m. . , where m. . is the ( i / j ) element of the 
1 j=i ^ 


matrix M. 
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Then 


s i 


4 [ B 0i p i-X ( X + Pi-1 V*' Vi)' 

+ Pi B u + (1-X)(B 0 . p ± _ x + B 2 p i+x ) 

+ P i( B Oi P i-X V X ' h i-l> + B 2i P i+A * 2 < X - h i>) 
+ B 2i P i+1 l X + P<- a 


( X + P i+1 a i< X ' h i>). 


For h i satisfying (5.73), we have S i > 0, which implies 


s i ih i A i>0 


where 


A i " B 0i Pi-A + Pi B li + B 2i Pi+X > 0 


N-l - 

Also, l m^. S. = 1 

i=l K ' 1 x 


k = 1(1) N-l 


1 th 1 

where m. . is the (k,i) element of the matrix M . 

K, 1 


Therefore 


N-l ! 

. ^ m k,i “ min sT 
11 1 s i < m 1 


h? A. 
x 0 x 0 


for some i Q between 1 


From (5.74) we have 


N-l 


= E < l i T < h i) 

3 i=l 3 ' 1 1 


which implies 


N-l 


|e.| * l |T(h.)| 

J i=l 


(5.75) 


(5.76) 


(5.77) 


(5.78) 

and N-l 


(5.79) 


j = 1(1) N-l 


(5.80) 
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with the help of (5.34), (5.78) and (5.80), it follows that 

I e jJ = °(h?) 

i. e. 

|| E | = 0(h 3 ) 

where h = max h^. 

In addition, if o'. £ 1 (uniform mesh), (5.34) gives 

T(h.) = 0(h?) 

Therefore, from (5.80), we have 

I E | = 0(h 4 ) 

where h — max h^ , i.e. our methods reduces to a family of 
order methods for uniform mesh. 


5.6 NUMERICAL EXAMPLES : 


Example 5 . 1 (Doolan et. al. 1980) 


ey' ' = y + cos nx + 2n e cos 2nx 
y(o) = y(i) = o 


The exact solution is given by 


Y(x) 


. , .-*•) . „„2 


l. * 


cos nx 


Example 5.2 (Carrier 1970) 

ey" - (2 -x) 2 y -1 
y(-l) = y(D - 0 


(5.81) 

(5.82) 

(5.83) 

(5.84) 
fourth 


We have taken Carrier's approximate solution for comparison purpose 
which is given by 
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with the help of (5.34), (5.78) and (5.80), it follows that 

l®il = °(h?) 

i.e. 

I E I = 0(h 3 ) 

where h = max h. . 

i 

In addition, if cr. s l (uniform mesh), (5.34) gives 

T(h.) = 0(hf) 

Therefore, from (5.80), we have 

« E 1 = 0(h 4 ) 

where h = max h. , i.e. our methods reduces to a family of 
i 

order methods for uniform mesh. 


5.6 NUMERICAL EXAMPLES : 

Example 5 . 1 (Doolan et. al. 1980) 

2 2 

ey' ' = y + cos 7TX + 2rr e cos 2 tix 
y (0) = y (l) = 0 


The exact solution is given by 


y(x) 


' -(l+x)/v'e + g-x/v'e] 2 

* ■■ 1 - ■ - 1 ■' -i — ■ - — — — 1 “ COS 




COS 7TX 


Example 5.2 (Carrier 1970) 

cy" = (2 -x) 2 y -1 

y (-1) - y(i) = o 


(5.81) 

(5.82) 

5.83) 

5.84) 
fourth 


We have taken Carrier's approximate solution for comparison purpose 
which is given by 
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y(x) = ( 2 - x 2 )~ 1 - exp [-(x+l ) /v'e ] - exp [-(l-x) /v'e] 

Example 5.3 (Carrier 1970) 

ey" =1 -2b(l-x 2 )y - y 2 

y (“i) = y(i) * o 

This example is solved for b = 0.0 and 1.0 by Ascher and Wiess 
(1984), for b = 0.5 by Roberts (1986). We have tabulated results 
for b = 0.0, 0.5 and 1.0 and taken Carrier's approximate solution 
for comparison purpose, which is given by 

y(x) = y R (x) + 12 exp(p x )/ [ l + exp (pi) ] 2 

+ 12 exp(p 2 )/ [ 1 + exp(p2) ] 2 

where y„(x) is the solution of the reduced problem given by 
y R (x) = -b(l-x 2 ) - [b 2 (l-x 2 ) 2 +1 ] 1/2 

and 

Pl , = (2/e) 1/2 ( 1 ± x) + 2 log(V2 + /3) 

Using the quasi-linearization method, we have 

c U-! = [- 2b(l-x 2 ) - 2 u n )u n+1 + U 2 + 1 

u n + l< 0) " W 15 = 0 
The initial guess u Q (x) is taken as 

U 0 (0) = u (1) = 0 and u Q (x) = y R (x) elsewhere and stopping 

criteria is 

max | u n+1 (i) " u n (i) I " ^ V 1 


where £ is a given constant. 
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5.7 DISCUSSION : 

We have implemented our method on three examples. Maximum 

errors at the nodal points, i.e. max | y. - y(x.)| are tabulated in 

i 1 

Tables 5. 1-5. 3 for different values of e, N, X and a. 

The second order method developed for general linear 
singularly perturbed boundary value problems using continuity 
condition can be applied here, but it will only give second order 
accuracy. Comparison between the results of previous chapter and 
this chapter for example 1 shows that error is reduced because of 
high order approximations. 
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Table 5.1 


Absolute maximum error in example 5 . l 








X 






N 

cr 

0.85 


0.90 


1.00 




1.00 

8.44233 

E -2 

2.62515 

E -2 

5.75677 

E -2 



60 

1.05 

2.27848 

E -3 

2.53310 

E -3 

1.08140 

E -2 

e* — I n 

-5 


1.15 

9.50132 

E -5 

9.98950 

E -5 

1.14810 

E -4 




1.00 

2.26680 

E -3 

4.90570 

E -3 

1.74408 

E -2 



120 

1.05 

3.75959 

E -5 

1.69450 

E -5 

1.10074 

E -5 




1.15 

2.99077 

E -5 

8.12125 

E -5 

2.15746 

E -5 




1.00 

1.34668 

E +1 

3.52120 


1.00872 




100 

1.05 

7.76747 

E -1 

6.63801 

E -1 

9.08492 

E -2 

/-.—i n 

-8 


1.15 

3.40082 

E -4 

3.31122 

E -4 

3.35928 

E -4 

o — XU 



1.00 

1.09771 


1.29771 


1.00427 

E -1 



200 

1.05 

8.35841 

E -4 

1.46160 

E -3 

5.47756 

E -3 




1.15 

2.47516 

E -5 

1.24349 

E -5 

1.51054 

E -5 




1.00 

3.43029 

E +1 

1.34162 


1.01017 

E -1 



150 

1.05 

1.12917 

E +1 

8.63424 

E -1 

9.30022 

E -2 


- 10 _ 


1.15 

4.20789 

E -5 

5.57506 

E -5 

5.46106 

E -5 

n 

II 

O 


1.00 

1.00075 


3.87275 


1.01007 

E -1 



300 

1.05 

4.85113 

E -4 

8.80132 

E -4 

3.28161 

E -3 




1.15 

2.20034 

E -5 

5.55807 

E -5 

1.49807 

E -5 




1.00 

2.23670 


2.15699 


1.01020 

E -1 



200 

1.05 

6.87588 


1.21240 


9.32885 

E -2 


-12 


1.15 

3.03721 

E -5 

1.99388 

E -5 

2.27321 

E -5 

C— 10 



1.00 

1.67117 


1.74504 


1.01020 

E -1 



400 

1.05 

3.21882 

E -4 

4.89111 

E -4 

1.92046 

E -3 




1.15 

3.02942 

E -5 

1.98478 

E -5 

1.50451 

E -5 
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Table 5.2 


Absolute maximum error in example 5.2 







A 





N 

cr 

0.85 


0.90 


1.00 



1.00 

8.27090 

E -2 

2.53335 

E -2 

5.61221 

E -2 


120 

1.05 

1.28270 

E -3 

1.38051 

E -3 

1.55133 

E -3 

„ i rt ^ . 


1.15 

1.37861 

E -3 

1.37754 

E -3 

1.37500 

E -3 

C— -LU 


1.00 

1.16652 

E -3 

5.90189 

E -3 

1.82053 

E -2 


240 

1.05 

1.38983 

E -3 

1.38973 

E -3 

1.38970 

E -3 



1.15 

1.38116 

E -3 

1.38026 

E -3 

1.37815 

E -3 



1.00 

1.02338 


2.16152 


9.89045 

E -2 


200 

1.05 

2.19239 

E -2 

8.88709 

E -4 

3.26500 

E -2 

i 

CO 

1 

o 

II 

(0 


1.15 

3.48469 

E -5 

3.39509 

E -5 

3.20879 

E -5 



1.00 

1.09756 


9.29963 

E -1 

9.94366 

E -2 


400 

1.05 

4.18056 

E -5 

4.17778 

E -5 

4.17120 

E -5 



1.15 

3.49235 

E -5 

3.40341 

E -5 

3.22187 

E -5 



1.00 

8.68514 

E -1 

9.70422 


9.96923 

E -2 


300 

1.05 

1.14491 

E -2 

2.22162 

E -3 

2.47163 

E -2 

- in” 10 


1.15 

8.87946 

E -6 

9.72226 

E -6 

1.22260 

E ~5 

e— io 


1.00 

9.19354 

E -1 

9.85432 

E -1 

1.00339 

E -1 


600 

1.05 

4.08246 

E -6 

4.06773 

E -6 

4.03352 

E -6 



1.15 

8.87869 

E -6 

9.72143 

E -6 

1.22250 

E -5 



1.00 

1.63941 


1.62333 


1.00022 

E -1 


400 

1.05 

5.94386 

E -3 

3.02194 

E -3 

1.81159 

E -2 

O in - 12 


1.15 

1.10274 

E -5 

1.20722 

E -5 

1.47466 

E -5 

e =10 


1.00 

3.46654 


2.76668 


1.00518 

E -1 


800 

1.05 

1.10274 

E -5 

1.20721 

E -5 

1.47466 

E -5 



1.15 

3.23433 

E -7 

6.41523 

E -6 

2.85756 

E -7 


Table 5.3 (a) 

Absolute maximum error in example 5.3 (b = 0) 
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X 



N 

O' 

0.85 


0.90 


1.00 



1.05 

1.61162 

E -3 

1.96165 

E -3 

2.74121 

E -3 


120 

1.10 

7.18562 

E -6 

7.49952 

E -6 

8.57211 

E -6 

1A ^ . 


1.15 

1.22547 

E -5 

1.30979 

E -5 

1.54331 

E -5 

e— xu 


1.05 

4.12840 

E -7 

4.32614 

E -7 

4.98784 

E -7 


240 

1.10 

2.43147 

E -6 

2.60394 

E -6 

3.10189 

E -6 



1.15 

1.12607 

E -5 

1.20405 

E -5 

1.43034 

E -5 



1.05 

5.63527 

E -3 

7.80762 

E -3 

1.12728 

E -2 


200 

1.10 

5.29125 

E -6 

5.57571 

E -6 

6.45412 

E -6 

i 

00 

1 

o 

rH 

II 

U 


1.15 

1.13972 

E -5 

1.22301 

E -5 

1.44895 

E -5 


1.05 

2.94392 

E -7 

3.11404 

E -7 

3.64386 

E -7 


400 

1.10 

2.41536 

E -6 

2.59594 

E -6 

3.09548 

E -6 



1.15 

1.12856 

E -5 

1.21109 

E -5 

1.43491 

E -5 



1.05 

6.88183 

E -4 

5.76028 

E -3 

4.78037 

E -3 


300 

1.10 

2.59192 

E -6 

2.78538 

E -6 

3.31420 

E -6 

o in" 10 


1.15 

1.12046 

E -5 

1.20324 

E -5 

1.43744 

E -5 

e— 10 


1.05 

1.61166 

E -7 

1.73394 

E -7 

2.07105 

E -7 


600 

1.10 

2.41532 

E -6 

2.60003 

E -6 

3.09477 

E -6 



1.15 

1.12035 

E -5 

1.20314 

E -5 

1.43732 

E -5 



1.05 

6.97031 

E -4 

2.18612 

E -3 

8.98200 

E -4 


400 

1.10 

1.15156 

E -5 

2.35158 

E -5 

3.10965 

E -6 

0 in -12 


1.15 

1.12824 

E -5 

3.06506 

E -5 

1.43658 

E -5 

e=10 


1.05 

5.49768 

E -6 

5.33567 

E “6 

1.97697 

E -7 


800 

1.10 

1.10224 

E -5 

2.95868 

E -5 

3.09112 

E -6 



1.15 

3.26594 

E -5 

3.35720 

E -5 

1.43658 

E -5 


120 


Table 5.3(b) 

Absolute maximum error in example 5.3 (b = 0.5) 


\ 


N 

< T 

0.85 


1.00 


1.05 

1.61062 E -3 

1.96065 E -3 

2.74020 E -3 

120 

1.10 

6.50927 E -6 

6.58689 E -6 

7.52364 E -6 


1 . 15 

1.14087 E -5 

1.22518 E -5 

1.45869 E -5 

-5 

o 

iH 

II 

to 

1.05 

6.63025 E -6 

9.32859 E -6 

6.50359 E -6 

240 

1.10 

6.50710 E -6 

6.57344 E -6 

6.52460 E -6 


1.15 

1.04508 E -5 

1.12305 E -5 

1.33829 E -5 


1.05 

5.63527 E -3 

7.80762 E -3 

1.12728 E -2 

200 

1.10 

5.29012 E -6 

5.57458 E -6 

6.45299 E -6 


1.15 

1.13962 E -5 

1.22291 E -5 

1.44886 E -5 

-n 

£=10 

1.05 

2.93323 E -7 

3.10335 E -7 

3.63255 E -7 

400 

1.10 

2.41442 E -6 

2.59489 E -6 

3.09443 E -6 


1.15 

1.12846 E -5 

1.21099 E -5 

1.43481 E -5 


1.05 

6.88182 E -4 

5.76028 E -3 

4.78037 E -3 

300 

1.10 

2.59191 E -6 

1.02274 E -5 

3.31419 E -6 

n 

1.15 

1. 12046 E -5 

1.20324 E -5 

1.43744 E -5 

e=10 10 

1.05 

1.61652 E -7 

1.73384 E -7 

2.07095 E -7 

600 

1.10 

2.41530 E -6 

9.29281 E ~6 

3.09476 E -6 


1.15 

1.12035 E -5 

1.20314 E -5 

1.43731 E -5 


1.05 

6.97031 E -4 

2.18612 E -3 

8.98200 E -4 

400 

1.10 

5.37324 E -6 

4.06362 E -5 

3.10966 E -6 


1.15 

1.13545 E -5 

1.21148 E -5 

1.43658 E -5 

£=10 12 

1.05 

9.87236 E -7 

2.75375 E -6 

1.97699 E -7 

800 

1.10 

2.42010 E -6 

2.16882 E -5 

3.09112 E -6 

1.15 

1.12824 E ~5 

3.89226 E -5 

1.43658 E -5 
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Table 5.3(c) 

Absolute maximum error in example 5.3 (b = l.o) 

X 



N 

0* 

0.85 


0.90 


1.00 



1.05 

1.59898 

E -3 

1.94900 

E -3 

2.72854 

E -3 


120 

1.10 

3.37277 

E -4 

1.21724 

E -5 

1.22124 

E -5 

n -10' 5 - 


1.15 

1.21174 

E -5 

1.21216 

E -5 

1.22741 

E -5 



1.05 

1.98066 

E -5 

2.95402 

E -5 

1.45904 

E -5 


240 

1.10 

5.71400 

E -5 

1.33572 

E -5 

1.30403 

E -5 



1.15 

1.21174 

E -5 

1.21218 

E -5 

1.22741 

E -5 



1.05 

5.63526 

E -3 

7.80763 

E -3 

1.12728 

E -2 


200 

1.10 

5.27643 

E -6 

5.56089 

E -6 

6.43930 

E -6 

o— io _ ^ - 


1.15 

1.13830 

E -5 

1.22158 

E -5 

1.44754 

E -5 



1.05 

2.79844 

E -7 

2.96857 

E -7 

3.49581 

E -7 


400 

1.10 

2.40139 

E -6 

2.58148 

E -6 

3.08101 

E -6 



1.15 

1.12714 

E -5 

1.20967 

E -5 

1.43349 

E -5 



1.05 

6.88183 

E -4 

5.76028 

E -3 

4.78037 

E -3 


300 

1.10 

2.59197 

E -6 

2.80660 

E -6 

3.31405 

E -6 

c=10~ 10 - 



1.15 

1.12044 

E -5 

1.20323 

E -5 

1.43742 

E -5 



1.05 

1.61522 

E -7 

6.20482 

E -7 

2.06962 

E -7 


600 

1.10 

2.41517 

E -6 

3.01709 

E -6 

3.09462 

E -6 



1.15 

1.12034 

E -5 

1.20313 

E -5 

1.43730 

E -5 



1.05 

6.97031 

E -4 

2.18612 

E -3 

8.98200 

E -4 


400 

1.10 

1.15156 

E -5 

2.35158 

E -5 

3.10965 

E -6 

p in" 12 


1.15 

1.12824 

E -5 

3.06506 

E -5 

1.43658 

E -5 

e-±u 


1.05 

5.49768 

E -6 

5.33567 

E -6 

1.97697 

E -7 


800 

1.10 

1.10224 

E -5 

2.95868 

E -5 

3.09112 

E -6 



1.15 

3.26594 

E -5 

3.35720 

E -5 

1.43658 

E -5 


CHAPTER VI 


Third Order Variable Mesh Methods For Non-linear 
Singularly Perturbed Boundary Value Problems 

6#1 introduction : 

Consider the following class of general non-linear singularly 
perturbed boundary value problems 

ey" = f(x,y,y') a < x < b, 0<e«l (6.1) 

y(a) = A, y (b) = B (6.2) 

where we assume that f is a smooth function satisfying: 

( i) — f (x,y,z) s o 

dZ 

(ii) — f(x,y,z) > o 
ay 

(iii) (— - — ) f(x,y,z) = oc > o 

ay dz 

and (iv) the growth condition 

f(x,y,z) = 0(|z| 2 ) as z — > » 

for all (x,y, z) e [a f b] x IR 2 

Under these conditions, the problem (6.1) -(6.2) has a unique 
solution, (Doolan et. al 1980) 

These problems occur in many areas of engineering and applied 
mathematics, notably fluid mechanics. Many practical areas where 
these types of problems occur and the difficulties associated m 
solving them have been discussed in previous chapt 
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In this chapter, the methods developed in the previous chapter 
for semi-linear problems are generalised for the non-linear problems 
of the above form. A two parameter family of third order variable 
mesh methods using cubic spline approximations is obtained, which 
gives more accurate third order approximations for first derivative 
terms also. These third order variable mesh methods reduce to a 
family of fourth order methods when the mesh ratio o\ is taken to be 
unity. Also when both o^and A are taken to be unity, a well known 
Numerov method is obtained as a particular case for the given 
problem with e = 1 and f independent of y' . 

6.2 DERIVATION OF THIRD ORDER METHODS : 

Let x Q - a , Xi = a +Yv i=l(l)N-l h* = V ^ = b. 

Given the values of y (x Q ) ,y (x^ , . . . ,y (x N ) of a function y(x) at the 
nodal points x Q ,x lf . . . ,x N , the interpolating cubic spline can be 
given in the following form (Ahlberg et.al. 1967) 


(x i+l" x)3 (x " x i ) 3 

s < x > = -SET M i + 6h . ~ M i + l + 


r h i H i 1 x i+l ' x 

L y( x i^ - -P ] < h. > 


h i 


+ 



h i M i + i 



X^ ^ X ^ X i+1 > i = 0(1)N-1 (6.3) 

where = S /7 ( x f) 
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and its derivative is given by 


j- (x - x i ) h i x. - x)“ h. , 

S ' (x) = [— ^7 6- ] M i + i + ^ ]«i 


“( * i+1 - h. 


y(x i+ i>- y(x.) 


x i a x - x i+1 . i = 0(1)N-1 (6.4) 


Now consider the following difference scheme for computing the 
approximate solution ' * * " ' ^N-l a t the nodal points 

X l /X 2 ' * " * ' X N-1* 

- e a oi y i-i + c a n y r e y i+ i + h l( B oi f i-x + B u v B 2i w = 0 


i=l(l)N-l 1//3 < X * 1 


(6.5) 


where and f i±A are approximations of f (x^,y (x i ) ,y' (x^) ) and 

f ( x i+X'Y ( X i+x) # y 7 ( X i+x) ) res P ec tively. The coefficient a's and B's 

are determined by replacing the approximate values by exact values 

lr (r) 

and equating the coefficients of h^y v ' (x^) , r = 0(1)4 to zero 
after taking Taylor's series expansion. 


Thus, as in previous chapter, we have 


a n . = cr . 
Oi l 


a n - 1 + 


( 6 . 6 ) 


(6.7) 


B 0i - 


o- 1 (<T i -l) ( 1-2X) +1 


12X 2 cr. 


( 6 . 8 ) 


B li = 


(cr?+ 1) (2A-1) + 2tr i (cr 1 +l)X(3X-l) 


12X 2 <r 2 


(6.9) 
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_ °i + (^-1) (2A -1) 

B2i HTV? 

1 


( 6 . 10 ) 


where cr^ = h^/h.^. 


i=l (l)N-l 


Consider the following approximations of the first derivatives at 

the nodal points x. , x*, ■ • ■ / x w _ in terms of v . v v 

* J 1' 1 2 ' ' y N“l* 

y i-l = h.(l+<r.) [ - y i+l + <°i + l) 2 Yi - ^.(cr. + 2) y._J 


i=l(l)N-l 


( 6 . 11 ) 


n+1 = hT liTo.) [ ( 2 <r i + 1 > y i + i " (°' i +l) 2 y i + ^Yi.! ] 


i=l(l)N-l 


( 6 . 12 ) 


and y i h i (i+<j i ) [ Y i+1 + ” °i ^i-i] " 6(o\+i)e [ f i+i“ f i-i] 


i=l (l)N-l 


Now we define f and f. as 
1 1±A. 


f i = f(x i' y i' y i> 


(6.13) 


(6.14) 


f i±A = f(X i±A' S(X i±A ),S ' (X i±A )) i=l(l)N-l (6.15) 


where S(x^ + ^) and S' (x. + ,) are obtained by putting 
x = x^Ah^ and x -Ah. . in S(x) and S' (x) , that are given by 


S(x) = 


(x i+r x > 3 .-. . < x - x i> 3 .-. 


M 1 + -6 h f - M i + 1 + 


r h i M i i , x i+i _ x , 

L y i - J ( KI 3 


‘i“l + l 1 x - x i 

6 J K h ± 


+ y. 


x i s x 3 x i+1 , i = 0(1) N-l 


(6.16) 
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(x - x,r h, 


[ \ A A -! ^ 1 r -(X. 

— JSI rKi * [ -§^ 


X) h_. 
+ 


y i+l" y i 
h i 


H s i 

x i" x " x i+l' 1 = 0 (1) N-l (6.17) 


where e M. = f( x j' y j ,y :p' 


1 = i# i+1 


After obtaining the values of f i# f i±x from Eqs. (6. 14) -(6. 15) and 

putting the values of a's from Eqs. (6 . 6) - (6.7) in Eq(6.5), we can 
write the finite difference scheme in the following form 


-e er.y.^t c(lHr ± ) y.- e Y i+1 + h.( B Qi f._ x + f l + B 0 , f 4A J = 0, 


2i i+X' 


(6.18) 


with 


y 0 = y(a) = A , y N = y(b) = b 

where B's are given by Eqs. (6 . 8) - (6. 10) . 

If the differential equation is linear, then the resulting 
tridiagonal system can be solved by Thomas algorithm, in the 
non-linear case, the system can be solved by quasi-linearization 
technique (Doolan et al.,1980) described in section 6.3. 


Remark 6 . 1 For c = X = <r^= 1 Vi, our finite difference scheme 

(6.18) reduces to well known Numerov's method for the problem 
y"=f(x,y), which is given as 

■ y i-i + y i+ i + + ‘ 10 f(x i' y i> + f(x i+i' y i+i ) ]"° 

(6.19) 
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6.3 MOW-LINEAR PROBLEMS : 

For non-linear problems of the form (6.1)-(6.2), we use the 
quasi-linearization technique. Starting from the initial guess 

u 0 “ u 0 (x) ' we obtain the sequence of linear boundary value 

problems 


eu 


n+1 


f y' ^ x,u n' u n^ u n+l + f y< X ' W u n+1 + f 
’ f y' ^ x,u n' u n^ u n " f y < x 'V u ;,> u n (6-20) 


u n+l (a) 


u n+l (b) = 13 


( 6 . 21 ) 


2 2 

It is known that (Doolan et.al. 1980), if 5_£( x 'y' z ) an( } fL_f(x,y,z) 

ay 2 ay 2 

are bounded for all x cs [a,b] and all real y and z.Then if 

a 2 f 

s o, there exists a constant C, independent of n and e, s.t. 

ay' 


u n+l (x) “ u(x) 


* C 


u n (x) - u (x) 


( 6 . 22 ) 


where u(x) is the solution of (6.1)-(6.2). 


6.4 ERROR ANALYSIS : 

Substituting for y^, the actual solution y(x^) in Eq (6.18), we 
have 

-e °‘ i y(x i _ 1 ) + e(l+cr i )y(x i ) - cy(x i ) 

+ h l [B oi f(x i . x ,y(x._ x ),y'(x i _ A )) + B U f(x 1 ,y(x.,y( Xl ) ( y'(x.)) 

+ B 2i f < x i+A'y< x i+A)'y' (x i+A ,) ] = T(h i } 1=1(1)N_1 


(6.23) 
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where is the local truncation error and defined by 

T(hi> = [ (5X " 3)(<r i " <7 i ) + 5X(2X " 1)(<r i ■ ^)] h fy lv ’ C x i> 


+ 0(h°) 


(6.24) 


Subtracting Eq(6.18) from Eq(6.23) and denoting e i = y(x i ) - y^ 
we have 


- E <r 1 e i _ 1 + e (l+<r.) ^ - e e ± 


h i { B 0! [ f ( x i-X'y( x i-X>'y'( x i-x)) - '(*i-x- S (*i-x)' 5 '( x x-x))] 

+ B ljL [f(x i ,y(x i ),y'(x 1 )) - f (x i ,y i ,yj L )l 

+ B 2i [^i+X'^W'^i+X^ - f ( x i + x-S( x i + x),S'( x i + x))]} 


= T(h ± ) 
Now let 


(6.25) 


y< x i±x> 

- b < x !±x> - (y< x i±x> - 

s(x i±x )) + (S(X i±X> “ § ( x i±X )) 



= e I± + e D± 

(say) 

(6.26) 

y' ( X i±x> 

- s'(x 1±x ) =(y'(x i±x )- 

s'(x i± x)) + (s'(x i±x )- s-(x i±x 

)) 


= e i± + e o± 


(6.27) 

and 

y' (Xj) - yj = e j 

j = i,i±i 

(6.28) 


From the theory of splines, it is known that(Ahlberg et. al. 19 67) 
e I+ = C + h? + 0(h B ) (6.29) 


e i± = C i h i + 0(h i> 


and 


(6.30) 
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where C ± are constants independent of h i and 

c; are constants independent of h i and has o^-l as a factor i.e., 

c; = D i(°i”l) where are constants independent of h.. 

Since, the approximations y}_ ±1 of y' (x i+1 ) are taken to be same as 
in chapter 4, (see eq. 4.32 & 4.34) we know that these are second 
order approximations, therefore we take 


e i-i 


e i+l 


where 


c i = 


c i h ? + °( h i> 

c 2 h l + °< h i> 


h? 

y'"(?3 ) + 


h? 


2 + <r. 


• r f £. T U . \ 

e-[( r+-?7 )y'"<’) 


CP. - 1 




(6.31) 

(6.32) 


and C 2 - 


2<r 2 

1 


<0 + 




2o* . +1 




For finding e^ = - y' (x^) , we do the following. 

Taking the usual Taylor series expansion for y around x^, we 
get the following expressions for y^ +1 and ^i-V 


i+1 


y i + h i y i 





IV 

y 



) 


(6.33) 
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h 2 

*i-i = *i - h i-i*i + r. 


h 3 

h i-i 


6~ y i" + sr 1 y iv (Hs n ) 

< i > , ( i ) 


x i-i<?5 €« ' <x i+1 

(6.34) 

Multiplying (6.34) by cr^, and subtracting it from (6.33), we get the 
following expression 

h 3 

*1+1 + ( <r i _1 >*i - cr i*i-i = h i( 1+<r i>yi + e? ( l+<r i)Yi" 

i 

h* 


3J ( - -ij y Lv (Ze‘’l ) 


(6.35) 


which implies 


* i+i + ^ r 1 '*! 


o-?v . 
x y i-l 


h id +cr i) 


*i + ^ y 


h i 

_ y 7 . 7 7 

6cr^ J l 


h ? 


24 (l+o* 


1+57) (** 


1V (Cg l) ) 


1 yiV ^ U ), 


Also, we have 


h? 


• • _ / 


l+l 


= y-'+ h i yi" + “I * (?s ) 


<r? 

1 


6 

(6.36) 

(6.37) 


h 2 

i-1 . iv/^d) 


and Y'i-! = yj/- h i_i y i" + -^- ± y ) 


(6.38) 


which implies, 


i_2 

h. , 


V 7 . 7 - \r' • — 

*1 + 1 


Y i-l “ V 1 + 57> *1 


y lv (^ u ) 


% y iv (€‘ u ) 

<T . 

1 


(6.39) 
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Therefore, 


tr.j+l)e [ f i+l f i-l] 6(0^+!) [ y i+i“ y j.li] 


-z-i-y'." + — A fv iv r? (i) \ 1 u iv /e (n,l 

6o i 1 12 (l+cr . ) [ y ^5 > ~2 y ^6 ) 

1 O' , J 

1 


Now subtracting this equation from (6.36), we have 

*1+1 + (<r l- 1)y i ~ ‘foi-i _ r 

h i( 1+cr i) 6(<r i +l)c |_ i+i~ f i-iJ 


(6.40) 


h/ 

y' - — i iv f f ( D \ 1 

Y 1 24(1+0-.) y l5 5 ' ~2 Y (? 6 M 

1 tr . * 


(6.41) 


which imlpies, 


h, ^ 

Y i y ^ X i^ 24(1+0'^) ( y ^5 ^ “ o .2 y ^6 (6.42) 

1 

which shows that y^ is third order approximation to y' (x ^) , 
therefore we can write 


e' = C 3 hl + 0(hf) 


(6.43) 


C 3 = 


(1+^) ( yl 


1V (^ U ) 


V <C> 

% 


where 
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From mean value theorem, we have 


£(x j ,y(x j ),y'(x j) )-f( Xj ,y j ,y' ) - §| x= (y ( x.)- yj) + g. 


- u^e^. + v^e^ , j=i,i±i 


(y # (^)-y^) 

X=Xj J J 

( 6 . 44 ) 


where 


af 

ay 


x=x. 


, af 
and ay' 


_ are denoted by u. and v. respectively. 

l — x • J j 


Subtracting ( 6 . 16 ) from ( 6 . 3 ), we have 
S(x) “ S(x) 


(x. +1 - x) h.(x. +1 - X) 


r 'n+i 
= [ 6hi 


] (M i " M i> + 


(x - Xl ) h i( x - X.) 


■ 

.-SET 


Substituting 


® ~~ ] (M i+r M i+i> + -h- 1 h x i+i> - ^i+i) 


( 6 . 45 ) 


and 


e M-, = f (Xj,y(X_.) ,y' (Xj ) ) 

E Mj = t (XyYj,y' .) , j = i, i+1 


and applying the mean value theorem, we have 


S(x) - S(x) 


(x.j.,- x) h (Xi - x) 


p A i+l 

- L 6h. 


1 ( } x i+l~ x 

J[ u i e i + V i e i) + T— "i 


6e 


(x - x i ) J h t (x - x t ) 


■ 

6h~ 

L l 


6e 




i+i e i+i + v i+i e i+i l + 


x - x. 


h. e i+l 


( 6 . 46 ) 
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Putting x = x i + Xh^, we get 
e - Sfx.+Xh^ - S(x.+Xh.) 


h? 


1 i 2 r 
— M* -1) u. 

6e L 


i+l e i+l + v i + l e ' 


ul 


h? 


i 2 r 1 

— (1-^) ((!-*) -1) ^ u i e i + + Xe i+1 + (6.47) 


Similarly putting x = x^-Xh.^ , we get 


'D- 


- S(x i -Xh i _ 1 ) - S(x i -Xh._ 1 ) 

m^ 2 -d + v i-i e i-i ] 


h2 i 

6ecr? 


h? 


((1-X) 2 -1) [ + v^l + Xe jL _ 1 + (1-X) ei 

P/T . L J 


6ecr^ 


(6.48) 


Subtracting (6.17) from (6.4), we have 


S' (x) - S' (x) 


r -(x. +1 - X) h. . y(x i+1 )- y 

" [ 2hT + S J (M i " M i> + hT 


r (x - x.) 2 h. -j 

[ 2hT 6 J (M i+1 M i+1* 


YCx^ - y l 


h i 


(6.49) 
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Again applying mean value theorem and substituting as above, we have 
S' (x) - S' (x) 


r -< x i + i- x) , h i 1 ( , 

= [ 2h.e 6e J ( u i e i + Vi 


r (X x L ) h x , , 

[ 2h i c 6 e J { u i 


i+l e i+l + v i+l e i+i ] 


e i+i" e i 
h i 


Putting x = x i + Ah^, we get 

h 


e i+l e i 


'D+ 


h 
h . 


+ ^ (A2 ' 1/3) [ u i + i e i +1 + v i+1 e i+1 ; 


^(1/3 - (l-X) 2 ) [ Ui e ± + v.ej 

Similarly putting x = we get 

e D- = S ' (x i~ Xh i-l> " S' (x i -Ah i _ 1 ) 

'i^i 


h. 

1 


h . 

i 

2 ccr 


(6.50) 


(6.51) 


fjH^-l/3) [ u.^e.^ + v i _ 1 e ^_ 1 


h. 

i 

2e<r 


i^(l/3 - (1-A) 2 )[u.e. + v.ej 


(6.52) 


Applying mean value theorem for the expressions in Eq(6.25), we 


have 
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•e + ed+o-.Je. -e e. +1 

i ( 8 »i[ 


+ h? 


at 

a y 


x=x i’ Ah i-i 


( y< x i-x> - sow) 


+ 


af 

dy' 


.( V ' |I W - S ' ( *i-x’) ] 


1 1-1 


+ Bli t ^Lx.l Y(Xi> ' y J + i'L x .( y ' (x i> - n) ] 


+ B 


r m 

2i L 3y 


x=x^+Ah^ 


( 


y(x i+x ) - i( X . +x ) 


+ 


af 

ay 7 


x=x^+Ah i 


( y'(x i+x ) - s'(x. +x) ) ] } = T ( h.) 


(6.53) 


Therefore, we have 

- c cr.e.^ + ofW.le, - e e . +1 

+ h i { B 0i[ U i-X ( e I- + e D-> + V i- A < 4- + e D-> ] 

+ B u [ u. e. + v. e; ] 

+ B 2i[ U i+* ( S I+ + e D+> + V i+X ( e I+ + e D+> ] } = T(h i> 


(6.54) 


where 



af 

i-A 

ay 


af 

i-A “ 

ay 


x=x r xh i-i 


' U l+A 


af 

ay 


x=X^+Atu 


* 


af 


x=x i- Ah i-i 


i+A dy' 


x=x i +Ah i 


Putting the values of e i+' e i+' e o+ and e^ ± and e^ (j-i,i±l) from 
(6,29)— (6.32) , (6.43), (6.47), (6.48), (6.51), (6.52), we get 
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- ctr i e i-i + E(1+0 i ) e i - c e i 


+ h i B oi 


h 2 


1+1 


h? 


( c - 1 * * it? ^ =t A ] 

L 1 J 


(l-X) ((1-X) 2 -1> [ u.e. + v.C 3 h 3 ] + x e._ x + (l-X)e. ) 


+ v i-x( c - h i 


, .3 ^ a i (e i- e i-l) h i 


h i 


^ (X 2 - 1/3) [u 1 _ 1 e i _ 1 + v._ lCl hf 


h. 

i 

2ccr 


y- (1/3 - (1-X) 2 )[ u.e, + v.C 3 h 2 ] ) 


+ B li [ U i e i + V i C 3 h i ] 


h? 


+ B 2i [ U i + x ( C + h t + “6c - X ^ 2 - 1 ) [ U i 


i+1 e i+l + v i+l C 2 hj 


n 


h i 2 r 3 i i 

+ gf- (1-X) ((1-X) 2 -l) [ U 1#1 + v.C.h 3 j + Xe i+1 (l-X)e. 


e_ . „ -e . h. 


f ' 3 c 1+ i e i “i 2 T 2 1 

+ v i + x l c + h l + h x + IF- -W [ “i + i e i + l + v i + l C 2 h i J 


h 

+ 2i < 3 / 3 - [ U i e i + V i C 3 h l ] 


+ O(h^) 


= T(h.) 


(6.55) 
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Clubbing the coefficients of e. e anH ^ ^ 

i-l' e i and e i +1 and shifting all 

otther terms to right hand side, we get 


H 'i + rir B 0i x(x2 " 1)u i-A u i-! + h i B oi u i-x X 


■ h i°i B oi v i-x “ TE5. B oi v i-x u i-i< x2 - 1 / 3 )] e i-i 


[ C(1+<ri) + B 0i (1-X * C ^ 1-x2 ) -1 ]u i « i _x + h 3 (l-X)B 0 .u... 


+ h B _ .cr. v. 
1 Ol 1 1 


i V i-X ■ 2^. B 0i V i-X (1 / 3 - (l-X)“)Ui + hjB llUl 


+ — — B 2i (l-A)[(l-l)-l,u.u i+x + ^(1-X)B 2 .u. +a - h.B 2 .v. 


'2C<T X B 2i V l+X^ 1 / 3 “ ( 1-A > > u i ] e i 


+ [' E + B 2 i A ( A2 -l)u i+ xU i+1 + h 3 B 2 .u i+A X + h.B 2 .v i+A 


+ -2F~ B 2i V i+X U i+l (x2- 1 / 3 >] e i 


= Tt^) - T 0 (h 1 ) 
= T(h.) (say) 


i=l (1) N-l 


( 6 . 56 ) 
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where T^I^) = 

u i-x ( C A + ~~2 MX 2 -1) V.^ c x 


h l 1 B 0i 


h? 


hf 


6ecr. 

1 


+ — ~ (1-X) ((l-A)-l) v. C.h? 

6 ecr^ 131 

h? 

+ v i-x( c - h i ■ ~ 2 ^r (* 2 - 1/3) v._ lCl 


h? 


2 h~ L - ( 1 -^) 2 V 3 ) ^ 


+ B. . v. C_h . 
11 l 3 l 


h? 


+ B 2i [ U i + X ( C A + SF- -V V i +1 C 2 


h 5 

+ g|- (i-x) ((i-x) 2 -i) Vi c 3 


h 3 


+ V. 

l+A 


/ o 11 J 

C h J + 

+ i 2e 


■d -1/3) + v i + l<= 2 


Clearly, 


+ 2? ^/ 3 * V i C 3 


T 0 (h.) = 0(h?) 


+ 0(h°) 


(6.57) 


Also, from (6.21) 


T(h.) = 0(h?) 


Therefore, 


(6.58) 


T(h j[ ) = 0(h?) 


(6.59) 
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putting Eq(6. 56) in the matrix-vector form, we have 

ME = T 


(6.60) 


where M is an (N-l)x(N-l) tridiagonal matrix whose ele m e n ts are 
given by 


m . . , ^ 

l, l+l 

= coefficient 

of e i + i 

in 

(6.56) 

i=l (1) N-2 

m . . 

1,1 

= coefficient 

of e . 

l 

in 

(6.56) 

i=l(l)N-l 

m . . , 

l, l-l 

= coefficient 

of e i-i 

in 

(6.56) 

i=2 (l)N-l 


tr 


E = (e 1 ,e 2 ,...,e N-1 ) ,T - (T (h ± ) ,T(h 2 ) , . . . ,5(1^ ) ) 


tr 


As in previous chapter, it can be seen that B B . and 

0i' li 

B 0 . are positive for 1/2 < X =s l if cr. satisfies 


(2A-1) 


1 

+ •— y 

/ 3 + 2X 

> 

V 

b 

V 

1 

+ 1 

/ 3 + 2X 

2 

2 V 

2X - 1 

d 

2 

. 

+ 2 y/ 

2X - 1 

d 


(6.61) 


Also it can been seen that , for 


h i < min(H 1 ,H 2 ,H 3 ) 


i=l(l)N-l 


(6.62) 


M is irreducibly diagonally dominant provided X € (1//3, 1]. 


where H 


_ u i-x u i-i ?l(x2 ~ 1) + 6 ecr i v i-i 
1 V 2c<r i X u._ x - v i . Jl u 1 _ 1 ) 


«2 = 


B 2i U i + A U i + l i < 1 -^> - 6EB 2i V i + i - 3B 2i V i-l U i + l (X ~ 1/3) 


B 2i XU i+l 
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H. 


[( 1_X) (B 0i U i-X + B 2i U i+x)’ B 2i v i+x3 6eo 'i +3 (l/3-(l-X) 2 )a-2B ,V 1 + .U. 

— — T — 1 21 1+A 1 


(1-X)[(1-X )-l](B oi U._ x u. - cr 2 B 2l u i+x u.) 


N-l 


Let S 


i- i -*■ 


i=l (1) N-l 


k=l 


(6.63) 


For hu satisfying (6.62), we have 


Si i 0 


and S i = h?G i + 0(h?) 


(6.64) 

(6.65) 


where G. = B oi u._ x + B^u. + B 2 .u. +X > 0 


Also, from the theory of matrices, we have 


N-i 


Y. m ji S i = 1 


j = 1(1) N-l 


i=l 


where m_^ is (},i) element of m” . 


( 6 . 66 ) 


With the help of Eqs. (6. 63) -(6. 65) and for ^satisfying (6.62), we 
have 


N-l 

X! “ji ' min S. “ ,2 
i=l i<i<N-l io io 


-1 < 1 
m . . - r- 


for some i Q between 1 & N-l 


(6.67) 


From (6.60) we have 


N-l 


e j " Y "S * (hi> 
1=1 


j = 1(1) N-l 


( 6 . 68 ) 



and therefore 
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i * h 

£ jl " Z2 


h iO G iO 


where h = max h. 
l=£isN 

Thus we have 


Ell - 0(h ) , 


j - 1(1) N-l (6.69) 


h = max h. 
lsi<N 1 


(6.70) 


Remark 6.2 If cr =1 Vi (uniform mesh ) , we have 
T(h i ) = 0(h®) , 

Hence we have 


||E|| = 0(h 4 ) , h = max h. (6.71) 

l^i^N 


6 . 5 NUMERICAL EXAMPLES : 

Example 6.1 : Convection - diffusion problem (Jain et. al.1984) 

e y" = y' 0 < x < 1 

y (o) = l, y(i) = o 


The exact solution is 


y(x) 


1 - exp(-(l-x) /e) 
1 - exp(-l/c) 


The boundary layer exists near the right end, so we take 
<r < 1, with h Q = (l-<r) / (l-o^) ,h^= ^i-1' where N is 

the number of mesh points in the interval. This ensures more mesh 
points in the boundary layer region. 



Example 2 : ( Kevorkian & Cole, 1981) 
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c y" = ( | - Dy' + § y o < x < l 

y(0) =0 , y (1) = l 

We have taken Nayfey's uniformly valid approximation, (Nayfey, 

1981) as exact solution for comparison, which is given as 

11 2 
y(x) = ~ 2 ex P(-(x - |)/c) 

The boundary layer exists near the left end, so we take <r > l, with 

N 

h^ — (o' “ 1) / (cr -1), h^= cr.jh^_^ , i=2,N .This ensures more mesh 

points near x = 0. 

Example 6.3 : (Bender & Orszag, 1978) 

ey" = -2y’ - e y 0 < x < 1 

y(0) = 0, y (1) = 1 


we have taken Bender & Orszag 's uniformly valid approximation as 
exact solution for comparison, which is given by 


y (x) = log - ( exp ( -2x/e ) ) log2 

The boundary layer region is at left end, so we follow the 
procedure as in example 2 for choosing cr. 

The initial guess is taken as u Q (x) = 0 and stopping criteria is 


u , . (x ) - u (x. ) I <10 
n+1 s i' n v i ; 1 


-7 


Vi 
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6. 5. DISCUSSION : 

We have implemented our method on three examples. Maximum error at 
the nodal points i.e., max |y. - y(x i ) | is tabulated in Tables 

6. 1 - 6 . 3 for different values of the parameters e, N, X, o*. 

It is observed that the errors are comparatively less for values 
of cr other than unity, which shows that variable mesh with 
approximations at off -nodal poiunts gives better results. 
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Table 6.1 

Maximum absolute error in Example 6.1 


X 



N 

<7 

0.85 


0.90 


1.00 



1.00 

9.81746 

E -1 

9.81746 

E -1 

9.81746 

E -1 


60 

0.90 

9.03708 

E -2 

9.25960 

E -2 

9.72461 

E -2 

C“10~^ - 


0.875 

5.78031 

E -4 

6.12345 

E -4 

6.84009 

E -4 



1.00 

9.80875 

E -1 

9.80875 

E -1 

9.80875 

E -1 


120 

0.90 

2.48247 

E -5 

2.66060 

E -5 

3.03009 

E -5 



0.875 

9.64444 

E -6 

1.03253 

E -5 

1.17449 

E -5 



1.00 

9.91992 

E -1 

9.91992 

E -1 

9.91992 

E -1 


150 

0.90 

9.99618 

E -1 

1.00000 


1.00004 


e —10 ® - 


0.875 

9.99826 

E -4 

1.05306 

E -3 

1.17058 

E -3 

b XU 


1.00 

9.99598 

E -1 

9.99598 

E -1 

9.99598 

E -1 


300 

0.90 

2.48001 

E -5 

2.65643 

E -5 

3.02596 

E -5 



0.875 

9.53883 

E -6 

1.02153 

E -5 

1.16196 

E -5 



1.00 

9.94998 

E -1 

9.94998 

E -1 

9.94998 

E -1 


200 

0.90 

6.92775 

E -4 

7.28057 

E -4 

8.01135 

E -4 

e=lo” 10 - 


0.875 

2.57257 

E -5 

2.71760 

E -5 

3.06257 

E -5 



1.00 

9.97499 

E -1 

9.97499 

E -1 

9.97499 

E -1 


400 

0.90 

9.50779 

E -5 

9.92227 

E -5 

1.10196 

E -*5 



0.875 

2.43328 

E -5 

2.72501 

E -5 

3.00868 

E -5 



1.00 

9.96363 

E -1 

9.96363 

E -1 

9.96363 

E -1 


300 

0.90 

7.80524 

E -4 

7.85154 

E -4 

7.86973 

E -4 

r in’ 12 


0.875 

2.09647 

E -4 

2.05907 

E -4 

2.12988 

E -4 

G— 1U 


1.00 

9.98181 

E -1 

9.98181 

E -1 

9.98181 

E -1 


600 

0.90 

7.77207 

E -4 

7.71633 

E -4 

7.71378 

E -4 



0.875 

4.42990 

E -4 

4.44399 

E -4 

4.93990 

E -4 
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Table 6.2 


Maximum absolute error in Example 6.2 







X 





N 

cr 

0.85 


0.90 


1.00 



1.00 

4.93354 

E -1 

4.93205 

E -1 

4.92889 

E -1 


60 

1.10 

4.92056 

E -1 

2.17676 

E -1 

4.80195 

E ~1 

-5 

P-10 - 


1.15 

3.33928 

E -2 

3.50604 

E -2 

4.91871 

E -2 



1.00 

4.92394 

E -1 

4.92176 

E -1 

4.91878 

E -1 


120 

1.10 

3.23356 

E -4 

3.25330 

E -4 

3.30103 

E -4 



1.15 

6.90069 

E -4 

6.92035 

E -4 

6.96799 

E -4 



1.00 

4.97903 

E -1 

4.97902 

E -1 

4.97901 

E -1 


150 

1.10 

4.99997 

E -1 

4.99997 

E -1 

4.99993 

E -1 

G — 10~® - 


1.15 

9.45442 

E -4 

9.73627 

E -4 

1.03259 

E -3 



1.00 

4.98946 

E -1 

4.98945 

E -1 

4.98942 

E -1 


300 

1.10 

3.20503 

E -4 

3.20507 

E -4 

3.20512 

E -4 



1.15 

6.87343 

E -4 

6.87342 

E -4 

6.87346 

E -4 



1.00 

4.98746 

E -1 

4.98746 

E -1 

4.98746 

E -1 


200 

1.10 

4.99999 

E -1 

4.99999 

E -1 

4.99999 

E -1 

-10 

e=10 


1.15 

3.44208 

E -2 

3.70527 

E -2 

3.21067 

E -2 



1.00 

4.99374 

E -1 

4.99374 

E -1 

4.99374 

E -1 


400 

1.10 

2.80828 

E -3 

2.86294 

E -3 

2.92001 

E -3 



1.15 

3.20223 

E -4 

3.20960 

E -4 

3.20548 

E -4 



1.00 

4.98997 

E -1 

4.99897 

E -1 

4.99897 

E -1 


300 

1.10 

4.99999 

E -1 

4.99999 

E -1 

4.99999 

E -1 

p in” 12 


1.15 

2.07622 

E -1 

1.62142 

E -1 

1.00568 

E -1 

C — 1U 


1.00 

4.99499 

E -1 

4.98997 

E -1 

4.99499 

E -1 


600 

1.10 

1.04976 

E -3 

1.07085 

E -3 

1.15266 

E -3 



1.15 

3.17444 

E -4 

2.85180 

E -4 

3.20599 

E -4 
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Table 6.3 


Maximum absolute error in Example 6.3 







X 





N 

O' 

0.85 


0.90 


1.00 



1.00 

6.75894 

E -1 

6.75894 

E -1 

6.75894 

E -1 


60 

1.10 

6.92318 

E -1 

6.92332 

E -1 

6.92324 

E -1 

-5 

p— 10 - 


1.15 

5.64459 

E -2 

5.65972 

E -2 

5.67671 

E -2 



1.00 

6.83299 

E -1 

6.82260 

E -1 

6.83234 

E -1 


120 

1.10 

5.90123 

E -2 

4.90086 

E -2 

4.58723 

E -2 



1.15 

4.02583 

E -2 

3.76094 

E -2 

4.10463 

E -4 



1.00 

6.84845 

E -1 

6.84845 

E -1 

6.84845 

E -1 


150 

1.10 

6.93147 

E -1 

6.34251 

E -1 

6.97642 

E -1 

C=10 -8 - 


1.15 

4.85552 

E -2 

4.26693 

E -2 

4.81394 

E -2 



1.00 

6.68986 

E -1 

6.68986 

E -1 

6.68986 

E -1 


300 

1.10 

9.10823 

E -2 

5.93021 

E -2 

4.75013 

E -2 



1.15 

4.27203 

E -2 

3.33572 

E -2 

4.18450 

E -2 



1.00 

6.88160 

E -1 

6.88160 

E -1 

6.88160 

E -1 


200 

1.10 

6.93144 

E -1 

6.93144 

E -1 

6.93144 

E -1 

e=lo" 10 - 


1.15 

4.38710 

E -2 

3.02597 

E -2 

4.18482 

E -2 



1.00 

6.90650 

E -1 

6.90650 

E -1 

6.90650 

E -1 


400 

1.10 

6.93114 

E -1 

6.93114 

E -1 

6.93114 

E -1 



1.15 

3.50238 

E -2 

3.23456 

E -2 

3.37818 

E -2 



1.00 

6.93256 

E -1 

6.93256 

E -1 

6.93256 

E -1 


300 

1.10 

6.84321 

E -1 

6.84321 

E -1 

6.84321 

E -1 

in -12 


1.15 

4.34522 

E -2 

4.24642 

E -2 

4.42711 

E -1 

c— ±u 











1.00 

6.93456 

E -1 

6.93456 

E -1 

6.93456 

E -1 


600 

1.10 

6.83243 

E -1 

6.83243 

E -1 

6.83243 

E -1 



1.15 

4.27654 

E -2 

4.42113 

E -2 

4.39877 

E -2 
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